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FOREWORD 


This technical report was written in the Applied Mathematics Section, Science 
and Mathematics Research Group, Warfare Analysis Department. Most of the work 
and all of the machine computations were done in 1974 and 1975. However, the 
essential mathematical steps of the present method were developed in 1968 shortly 
after the publication by D.R. Weidman and M.P. Brunnor of their model for the 
solution of the same problem. See NWL Technical Note TN-K/21-68. n. } ' ’ 
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ABSTRACT 


: ae A solution by deterministic methods is described of the problem of computing 
q the single-shot kill probability of a point target at a random point from a uniform 
A ; distribution over the interior of an arbitrary ellipse in the plane, given that the distri- 
bution of shots is uncorrelated bivariate normal with respect to a rectangular coordi- 

nate system in the plane, and that the weapon has a cookie-cutter damage function 

with prescribed lethal radius R. This solution has been programmed at NSWC, 

Dahlgren Laboratory. The numerical evaluation of a double integral, whose integrand 

contains the so-called elliptic coverage function, is required. Computer results clearly 

show the superiority of this solution over a non-deterministic, Monte Carlo method 

of Weidman and Brunner. 
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1. INTRODUCTION 


The problem treated in this report is as follows. A bomb falls in an uncorre- 
lated bivariate normal distribution, with standard deviations ¢, and o, and mean 
point at the origin, with respect to a rectangular Cartesian coordinate system Chk. 
The bomb is assumed to have a known lethal radius R and a cookie-cutter lethality 
function; that is, everything is destroyed in a circle of radius R with the impact 
point as center, and nothing is damaged outside this lethal circle. This lethality or 
damage function has been used in many weapons effectiveness studies. 


An ellipse of semi-axes a and b, a > b, is given, with center at a known point 
(u, v) in the Chk system and the major axis making a known angle ¢ (phi) with the 
h-axis (see Figure 1). A point target is at an unknown random position in the 
interior of the ellipse, governed by a uniform distribution over the area of the 
ellipse. It is required to find the probability P, that the point target is destroyed by 
a single bomb. 


An expression for P,, is derived in terms of a double integral, (2-12) or (2-13), 
which must be evaluated numerically. The expression is a function of eight variables, 
u, v, @, a, b, R, o,, 0,, of which seven are independent. We could normalize with 
respect to R, or o,,, or in other ways, and thus reduce the number of variables by one. 


The variables of integration in (2-12) and (2-13) are r and 0, parametric coordi- 
nates in the ellipse. In each element of area, the function P(h, k, CT) R) is 
evaluated by NSWC/DL subroutine PKILL, tho coordinates (h, k) of the impact 
point of the bomb being expressed in terms of u, v, ¢ (= b/a), ¢, rand @ (211), 
This subroutine is based on [2], and computes the probability that a single bomb, 
falling in an elliptical normal distribution of the type described above, will fall in a 
circle of center (h, k) and radius R. 


In special cases, namely when and only when o, =0,, we can use NSWC sub- 
routine CIRCV, computing the circular coverage function [5], instead of PKILL. 
CIRCV is at least ten times faster than PKILL. Subroutine CIRCV computes the 
probability P(R, d) that a single bomb, falling in a standard circular normal distri- 
bution, ie, o, =o, = 1, will fall in a circle of radius R with center offset a distance 
d from the mean point of the distribution. Further details on this special case are 
given in Appendix A. 
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This method, referred to hereafter as the DJH method, in cont: ust to a different 
approach discussed below, and in [1], is entirely deterministic, The accuracy of the 
computed P,, is limited only by that of the numerical evaluation of the double 
integral. If an iterative numerical double integration routine is used with successive 
halving of the step sizes in each direction, P, can be computed to any reasonable 
accuracy level by continuing the iterations until successive iterates agree with each 
other within a preassigned tolerance. A FORTRAN listing of tne program is given 
in Appendix F. 


A non-deterministic, Monte Carlo method is used in [1]. A simulated bomb 
impact point (h, k) is determined by the use of random normal numbers. One 
then computes the area common to the ellipse and the lethal circle, ie., the circle 
of center (h, k) and radius R. The conditional probability of destruction of the point 
target, given that the bomb hits at (h, k) is then given by the ratio of this common 
area to the total area of the ellipse. This follows from the uniform distribution of 
the target positions throughout the interior of the ellipse. By replicating this proce- 
dure a large number of times and taking the mean of the conditional probabilities, 
an approximation p is obtained to the overall kill probability P,. 


The method used in [1] for estimating the error in computing P,. is to use the 
probable error (P.E.). However, their expression for this quantity is not correct. The 
correct expression is given by 


pz. = 0.6745. / 2 =P) (1-1) 
n(n — 1) 


where p, Gi = 1,2,....n) is a conditional probability as described above, and f is the 
mean of these n conditional probabilities. This mean, f, is the estimate obtained 
for the kill probability, P,, by the particular simulation involving the n conditional 
probabilities. 


The validity of (1-1) is established in Appendix C. Although this corrected 
formula for the P.E., being smaller than the incorrect formula by a factor I Jn — | 
(see (C6)), puts their method in a better light, the computer results in Section 3 
show, not surprisingly, the superiority of the method (DJH) of this repost for 
computing P, . 


See Appendix C for a detailed discussion of the P.E. in the Weidman-Brunner 
(WB) method. 


pe ay Oa Ea et Se a a ere et ell ed do 


2. DESCRIPTION OF DJH METHOD 


The geometry of the problem is shown in Figure 1. A bomb falls in an uncorre- 
lated bivariate normal distribution centered at C, the origin of the rectangular 
Cartesian coordinate system Chk, with standard deviations o, and oc, along the h- and 
k-axes, respectively. In gencralo, > o,. Ifo, =o,, we have a special case (Appen- 
dix A). The probability density function of the impact point of the bomb is 


f(h, k) = 


Figure |. Geometry of the Problem 


7 : It is assumed that the bomb has a cookie-cutter lethality function with known | 
zy lethal radius R (see first paragraph of Section 1). | 


A point target is located at a random point from a uniform distribution over 
the area of a given ellipse, centered at a point O with known coordinates (u, y) in 
the Chk system. The semi-axes are a, b, with a> b in general. If a = b, the ellipse 
is a circle, and this is a special case (Appendix A). The major axis of the ellipse is . 
inclined at a known angle ¢ to the h-axis, measured positive counterclockwise from 
the h-axis and in the range - 90 < ¢ < 90 degrees, 


A rectangular coordinate system Oxy is set up with origin at O, the center of 
al : the ellipse, and x- and y-axes along the major and minor axes, respectively. 


We shall need the equations of transformation from the Oxy to the Chk 
coordinate system. We set up an intermediate system Ox'y’ (Figure |), with the 
x’- and y’ayes purallel to the h- and k-axes, respectively. Then the equations of 
translation between the Chk and Ox’y’ systems, and of rotation through the given 
angle @ from the Ox'y‘ to the Oxy system, are the following: 


f , h 


=x’ +u 
nr : pegs y (2-1) 
i sinners 
y’ = x sin @ + y cos ¢$ (2-2) 


Combining (2-1) with (2-2), we get the following equations of transformation from 
the Oxy to the Chk system, in terms of the constants u, v, and ¢. 


ee) 
k=vtxsing +ycos¢ (2:3) 


The equation of the ellipse (Figure 1) in the Oxy system is 


2 2 
cae (2-4) 


Let b/a = c, with O< cc < I except in the special case In which the ellipse is a circle 
(Appendix A), when c = 1. 


NG 


We represent points on or inside the ellipse by curvilinear coordinates, r, 0, 


i where 

f 

i x =r cos 0 | 

ij. y = cr sin @ (2:5) 
i 
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The curvilinear coordinates, r, 8 can be given geometrical significance in the 
following way. Since b = ac, (2-4) is equivalent to 


y? 


x? Lies = a? (2-6) 
and if we put 

: : vi (2-7) 
we have 

Btn? =o? (2-8) 


or, the ellipse corresponds to a circle of radius a in the &7-plane; while for interior 
points in the ellipse, (2-5) and (2-7) give 


ead, 
n=rsin @ 


(2-9) 


or, r and @ are polar coordinates on and inside the auxiliary circle represented by 
(2-8). 


By a standard theorem in analysis (for example, see (3, p. 348]), the element 
of area dA or dx dy transforms in the following way 


ax ay| 
3 a a 
da = (OY ge ag e [PT OTT ae ae 
a(r,é) ax oy 
Te) 


cos @ c sin 


~rsind cr cos 6 dt d6 = or dr dé 


(2-10) 


\a(x,y)/8(r,0)| being the Jacobian determinant. 
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Subroutine PKILL, or the elliptical coverage routine, [2], does the following 
calculation, Suppose we are given a rectangular coordinate system Chk and a bomb 
falling in an uncorrelated bivariate normal distribution with mean point C and standard 
deviations o, and oc, as in the present problem. Suppose we are alse given a circle 
with the center at a given point (h, k) and a given radius R, This subroutine computes 
the probability P(h, k, 0,, 0,, R) that a single bomb falls in the given circle, This 
probability, computed by subroutine PKILL, is also designated as P,,,,, but Pyiry 
should not be confused with the overall single shot kill probability P, computed by 
the DJH or the WB method of this report. Suppose that the lethal radius of the bomb, 
with the so-called cookie-cutter damage function, is also R, say 50 feet (this of course 
is the reason for taking this constant R as one of the input parameters in computing 
P(h, k, On» xs R)). Then a point target at (h, k) is destroyed if, and only if, the 
impact point of the bomb is in the given circle centered at (h, k). It is obvious from 
a simple geometrical figure that the point (h, k) is in the lethal circle of the bomb, 
if, and only if, the bomb’s impact point is in the circle of radius R centered at (h, k). 
Thus, the single shot probability of destruction of the point target at (h, k) is given 
by P(h, k, o,, o,, R) as computed by subroutine PKILL, 


We proceed to set up the double integral whose value, designated as P,, is 
the probability of destruction of the point target at a random point from a uniform 
distribution over the area of the given ellipse. Given any element of area cr dr dé 
of the ellipse (2-6), the probability that the point target lies in this element of area 
is cr dr d@/(rab), the ratio of the area of the element to tho entire area of the ellipse, 
since the distribution of the target position is uniform over the ellipse. If the target 
does lic in this element of area, and if a representative point in the element has 
nates, neglecting infinitesimals of higher order), then the conditional probability of 
destruction of the target is P(h, k, o,,0,,R)+ cr drdé/(mab). The coordinates of the 
point being (x, y) in the Oxy system, or (r, 0) in curvilinear coordinates, the (h, k) 
coordinates of the point, by (2-3) and (2-5) are given by 


h=utrcos@ Eamets . (211) 
k=vtrcosé@ sing ter sin @ cos¢ ‘ 

u, ¥, and ¢ being constants. Hence the conditional probability of destruction of the 

target, given that it Hes in the element of area cr dr dé, is | 


P(u +r cos 9 cos ¢ -- cr sin @ sin ¢, v +r cos 0 sin \ 
+ cr sin @ cos ¢,o,,0,, R) cr dr d@/(rab), 


which could be written more briefly as P(h, k, 0,» 0,)R) cr dr dé /(rab), with h, k 
given by (2-11). 


The overall kill probability P, is obtained by integrating these conditional 
probabilities over the entire ellipse, since the target must be somewhere in the 
ellipse. The curvilinear coordinate r goes from 0 to a (every point on the ellipse 
having this coordinate r = a), and @ from 0 to 27. There is no symmetry here, in 
general, because of the arbitrary position of the ellipse in the Chk system (Figure 
1), so that 9 niust be integrated through all four quadrants, 0 to 27. The result is 


1 dn pa 
PL = b ff Pou + cos 8 cos ¢ — or sin 0 sind, Vv trcos 6 sin} 


K na? 
+ cr sin 0 cos ¢,0,,0,, R)+ r dr dé (2-12) 


Here we have replaced c¢/(rab) by c¢/(za * ac) or (1/7a*). This can be written more 
compactly as 


| on ry * 
P, -— f fr PO borage RD dr do (2-13) 


with h, k given by (2-11). The P functions on the right are computed by subrou- 
tine PKILL. The double integral as written in (2-12) or (2-13) is now ready for 
evaluation by a numericz! dout’> integration subroutine, which will treat the variables 
rand @ as dummy variables like the variables of integration in any double inte- 
gral, so that their geometrical significance as curvilinear coordinates in the ellipse, or 
polar coordinates in an auxiliary circle, is no longer relevant to the calculation. 


One fact that we can deduce from (2-12) or (2-13), and which is of some 
theoretical significance, and is useful in designing test cases and criteria for deter- 
mining cases in which P, ~ 1 or Py ~ 0, is that the computed value of the kill 
probability, P,, is simply the space average, over the ellipse of integration, of the 
values of the elliptical coverage function P(h, k, 0,» 0,» R), appearing in the integrand 
in (2-13) and computed by subroutine PKILL, [2]. 


*See note, page 29, regarding h and k. 
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We can prove this in the following way. Let F denote the average value of 
P(h, k, o,, 0,, R). Then 


Ps Sf fr k, o,, 9,, R) aa f fan (2-14) 


where each double integration, in the parametric variables r and 6, is over the ellipse 
of integration. Since dA = cr dr dé by (2-10), we then have 


: J fre k, 0, 0,, R) cr dr do 
p= : oe (2-18) 
Sf dr dé . 


Now we move the two factors of c outside their respective integral signs and cancel 
them, and trivially evatuate the integral : 


2n s*a 
f | 1drd@ as 7a?, 
0 (¢) 


so tiiat we have 


1 an fa 
Be. ff Jf PO Ks oy 0) RD Fd a (2-16) 


which is the same as the expression for the kill probability P, in (2-13), so that we 
have 


(2-17) 


as stated. 


3. NUMERICAL RESULTS 


The basic equation for the computation of the kill probability P, by the DJH 
method is (2-12), or equivalently, (2-13), with h,k given by (2-11). We have to eval- 
uate numerically a double integral, and we have to call the subroutine PKILL, com- 
puting the elliptical coverage function P(h, k, 0 n? 9p R) [2], for each point used in 
the integration. 


The double integration routine used is a CDC software routine, using Simpson’s 
rule in each dimension, provided by Control Data Corporation, manufacturer of the 
CDC 6700 machine, This routine allows for the selection of different tolerances for 
the horizontal and vertical integrations. However, the authors of the present report 
did not give systematic study to the selection of these tolerances until several of the 
cases in the testing program required an excessive amount of computing time, more 
than one minute in some cases (page 16). The error analysis which was then worked 
out, as well as a detailed description of the integration method, is given later in this 
section, pages 21-29. Hore we give a brief general desoription of the method, assuming 
that the positive tolerances, denoted as ¢, and ¢,, have been appropriately chosen. 


We follow here the notation of (3-1) and Figure 3 (page 22), the integrand 
F(x, y) being unspecified und the integral being regarded as an iterated one. For each 
integration point y = s used in the vertical integration (with respect to y), the single 


integral fru, s) dx or f(s) is approximated by Simpson's rule, 


We give first an overall summary of the two integrations, followed by a few more 
details, but for fuller details ne pages 21-29, For a given y = s, where cg s< d, the 
exact value of the integral I F(x, s) dx is denoted by i and the Simpson approxima- 
tion by I. The calculation of I, involves a tolerance which is denoted here ase,. A 
function of y, namely I, or l,, is thereby defined on [c,d]. This function, in general, 
has discontinuities at values of y for which the number of subintervals of {a, b] 
changes, say from 16 to 32, subintervals (or from 17 to 33 points), but because of the 
smooth and slowly changing integrand here, (2-12), we can ate es the points of 
discontinuity are finite in number and, therefore, the integral VB l, dy is approxi- 
mated by Simpson sums [(d ~ ¢)/3n] & k, I, as the number of integration points be- 


comes infinite, the multipliers k, being the familiar Simpson factors 1,4,2,4,...,2,4,1. 
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An overall tolerance ¢, sey ¢ = 0.005, is specified. The objective is to use suffi- 
cient points in both the horizontal and the vertical integrations so that we have a high 
degree of confidence, on heuristic grounds at least, that the final approximation I is 
within ¢ of the true value of the double integral; or, |I — Ti< ¢. The tolerance Ey 


referred to above, and a tolerance e , used in the vertical integration (with respect to | 


y), depend on ¢ in a manner indicated below. 


For a given y, or s, sufficient points on [a, b) are taken to give a high degree of 
confidence that the approximation I, is within e, of the exact value I, or J DECK, s) dx. 


The usual rule in practice is that two successive Simpson sums aa z k, F(x,, 8), 
and a z k 5 F(x, 8) (where n, = 2n,, the x,’s, i= 0, 1,. veo, a nai spaced 


on i ; i the x, ’s,j ™ 0, 1,...,n, also equally spaced at half the mesh length of the 
x,'s) agree with each other within e,. A stricter criterion could be uscd if desired, for 
example, that of three successive Simpson sums, the largest and smallest must be 
within e, of each other. But, once a criterion has been selected, it will be applicable 
for every value of y ors on [c, d), and a function I, is defined on [c, d) as described 
above, Note that different numbers of points on {u, »), in general, are used on differ- 
ent horizontal lines, i.e., for different values of y or s, This is already implied by the 
remarks above regarding discontinuities in I, at values of y for which the number of 
points on [a, b] changes. 


It is shown that the integral fe I, dy differs by less than (d — c) é, from the 


exact value of the double integral, rh dy or T (3-8). But a sufficient number of 
points on [c, d] are taken so that (to a high degree of confidence at least) the approxi- 


mating Simpson sum sae xk, I, , oF das k, I, ) is within e, of the value of the 


integral ah l, dy (3-10), where e,, like €y» is one of the ar depending on the 


overall tolerance e. Therefore, the final Simpson approximation, 2 an Ly k, Le in the 


vertical integration, is within (d — c) €, +¢, of the exact value T of the double 
integral (3-13), 


This sum, (d — c) €, te,, is equated with the overall tolerance ¢ (3-15). The 
tolerances €, and e, may be chosen in any manner such that this equality is satisfied. 
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The most obvious choice is e, = €,, in which case each of these tolerances has the 
value e/[1 + (d -- c)] (3-17). 


We now consider in slightly more detail the computation by Simpson’s rule of 
a typical one of the horizontal integrals (i.e., with respect to x) f F(x, s) dx, where 


s is some fixed value of y on [c, d]. If n subintervals of [a, b] are used, or n+ 1 
points, n being an even irteger, the Simpson's rule approximation is 


b—a 


h 


n ba 
em 2K, F(X, 8) Fo 8) + 4F(a + ,$) 
3n '* 3n 


+ 2F(a + 2 aes s)teo t+ 4Fla + (n -- eam 8] + F(b, >} 
n n 


The step size (b — a)/n is repeatedly reduced by successive doubling of n until we reach 
a value such that Mo Pees Isnt < ey) in case the usual criterion described above is used. 
A similar method is used in the vertical integration, with tolerance e,. 


First we give the numerical values taken by the input parameters 0,» 0%)» R, u, 
v, a, b, and ¢ in the Table E3 in final form as they appear in this technical report, 
pages 61 through 111. Actually, values are assigned to the ratio b/a rather than to b 
itself, Later in this section, we discuss the testing program involving both the DJH 
and WB methods, and such topics as simple tests for detecting cases in which P, a 0 
or { without performing the double integration of (2-12) or (2-13), and the error 
analysis for the method of double integration as used in the DJH program, in the 
general case where P, is not near Q or 1. 


The standard deviation o, takes the value 1 in every case because of the normali- 
zation which is used, a takes the velues 0.1, 0.5, 1, and 2 on every page; the ratio b/a 
always takes the values 0.2, 0.5, 0.8, and 1; and the angle ¢ takes the values 0°, 45°, 
90°, and —45° on every page. Therefore, the only parameters whose values can vary 
from page to page are R, o,, u, and v (see Table 1), and these are the parameters 
whose values are listed in the appended table. The total nuinber of cases in the main 
table (Table E3) is 59,904, on 51 pages, 61 to 111. 


11 


' Table 1. Values of R, o,, u, and v for Main Table (E3) 


Values 
of 


Values of u Values of v 


: 0,2,0.5,0.8,1 0,0.5,1,1.5,2 0,0.25,0.50,1 - | 

q 0.2,0.5,0.8,! 0,0.5,1,1.5,2 0,0.25,0.50,1 | 

g 0.2,0.5,0.8,1 0,0.5,1,1.5,2 0,0.25,0.50,1 | 
| 0,2,0.5,0.8,1 0,0.5,1,1.5,2,2.5,3,4,5 0,0.25,0.50,1 
0.2,0.5,0.8,1 0,0.5,1,1.5,2,2.5,3,4,5 0,0.25,0.50,1 

0.2,0.5,0.8,1 0,0.5,1,1.5,2,2.5,3,4,5 0,0.25,0.50,1 

(See Note) | 

0.8,1 0,0.5,1,1.5,2,2.5,3,4,5 0,0.25,0.50,1 : 
0,2,0.5,0.8,1 0,0.5,1,1.5,2,2.5,3,4,5 0,0.25,0.50,1 

1 0,0.5,1,2,2.5,3,3.5,4, 0,0.25,0.50,1 ; 


4.5,5,5.5,6 


Note: For R = 2, 0, = 0.2: v = 0, 0.25, 0.50, 1.50 
But for R = 2, 0, = 0.5, 0.8, 1: v = 0, 0.25, 0.50, 1 


The tolerance, or ¢ in the error analysis, is generally 0.005, indicating that, 
although the main table gives values of P, to three decimal digits, the third digit { 
may be incorrect by several units, and we should round off to two decimal digits to 
obtain reliable results. The tolerance was taken as 0.001 in two pages of the table 
(pages 61 and 65) where o, takes its smallest value, 0.2, and R its two smallest 
values, 0.25 and 0.50. But a tolerance of 0.001 entailed the use of too much compu- 
ter time in general. 


At the top of each page are the vulues ofc, , o, and R for all of the cases on that s 
page, but since o, ™ 1 always, the significant values are those of o, and R. The values " 
of u, v and ¢ appear at or near the left end of the lines, and those of a and b as column 
headings. 
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The angle ¢, on every page, takes the four values listed above (--45°, 0°, 45°, 
90°), and similarly a and the ratio b/a take the four values for each which are listed, 
The most common values for u and v are u = 0, 0.5, 1, 1.5, 2, and v =0, 0.25, 0.5, 1. 
The larger values of u (greater than 2.0) are used with the larger values of R, since 
for these values of R and small to moderate values of u, the values of P, are generally 
near 1, The kill probubility P, increases strictly monotonically with R for fixed 
values of the other input parameters, and we always have P, + 1 as R*. On 
page 109 of the table, where o, = 1, R= 4, the smallest value of P, foru< 1 is 0.941, 
and for us 2 it is 0.838, the smallest value on the page. Hence, in order to cover a 
substantial part of the range (0,1) for P,,, an additional page (page 111) was computed 
foro, = 1, R = 4, with u taking values up to 6, 


Also, on pages 89-90 of the tables, the values of v are 0, 0.25, 0.5, 1.5, since on 
these pages, where o, = 0,2, R = 2, the values of P, are generally insensitive to the 
value of v, and we get a wider range of values of P, with these values of v than with 
the usual values 0, 0.25, 0.5, 1. 


Further comments on the selection of values for the input parameters are given 
later in this section, pages 18-20, 


In checking out the DJH method of this report (see (2-12) and (2+13)), 22 cases 
were run on the CDC 6700 computer. These cases are listed in Table 2. Actually, 
one case was inadvertently assigned two case numbers, 14 and 15, so that there were 
only 21 different cases, but we foliow the original numbering here. Some of these 
cases were checked by three different methods, (a) double integration using subrou- 
tine PKILI. (ie., the general DJH method), (b) double integration using subroutine 
CIRCV instead of PKILL (see Appendix A, puge 31), and (c) single integration using 
CIRCV (see (A2) and (A3) in Appendix A), Methods (b) and (c) can be used if, and 
only if, 0, = o,, as stated in Appendix A. 


The results were uniformly consistent, justifying the conclusion that the DJH 
program computes P, with an absolute error not exceeding a few units in the third 
decimal digit. Although this level of accuracy is adequate for most applications, five 
to six decimal digit accuracy can be attained for P, by a simple change in the PKILL 
subroutine to yield six decimal digit accuracy instead of three as it now stands, and a 
willingness to pay the costs of significantly increased computer time. The time will, 
in general, be more than doubled, in obtaining results correct to six rather than three 
decimal digits. Subroutine PKILL, [2], on the CDC 6700 machine, requires approxi- 
mately 5 milliseconds for three-digit accuracy and 10 milliseconds for six digits. In 


Table 2. Test Cases for DJH and WB Methods 
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addition, there will be a finer tolerance, say « = 107°, in the error analysis later in this 
section, pages 21-29, in computing to six digits, than the toler.nce of 107? in the 
three-digit calculation, and this will entail more calls to the PKILL subroutine. 


In case 16, computed by the three methods mentioned above, with the objective 
of obtaining six-digit agreement, the smallest and largest of the three results for Py 
were 0.02251899 and 0.02251922 respectively. We note that method (c), ‘the third 
method used, js quite different from the other two in setting up the integration. 


All except cases 8 through 12 were run by the method of [1], which will be 
referred to as the WB method (Weidman and Brunner). These five cases were designed 
to test the DJH method by using a sequence of decreasing ellipses closing down on 
the point (0,4), the smallest one having a major axis of 0.0625. 


The results obtained by ‘ice WB method were of variable quality. In some cases 
they were surprisingly good, as in Case 6 in which the DJH and WB methods gave 


" results of 0.019962 and 0.019422 respectively. An isolated case of good agreement 


could be due to chance, but in several other cases the agreement was quite good. 


But in some other cases the WB method, when applied more than once, generat- 
ing different simulated bomb impact points each time, gave results that differed signif- 
icantly from one another. Thus, in Case 2, which was run threc times with the WB 
and once with the deterministic DJH method, the WB method, using 25, 50, and 100 
simulated impact points, gave results of 0.011872, 0.005936 and 0.0217015 respec- 
tively. ‘Yet the apparently correct P,, as given by the deterministic DJH method, 
was 0.044321. Moreover, the 50% confidence interval for the WB mothod, computed 
by the meihod described in Appendix C (see also a brief description in the next para- 
staph of this section), in none of the three cases contained the correct value, 0.044321. 
For the calculation in which 100 impact points were used, the probable error was 
0.008998, and so the $0% confidence interval was 0.0217015 + 0.008998 or 
(0.012704, 0.030700). 


The 50% confidence interval is (B — P.E., B + P.E.), where Dis the mean of the n 
conditional probabilities { P| in a particular simulation, and P.E, is the probable error 
computed by (1-1) or (C5). For a justification of this expression for the 50% confi- 
dence interval, see the latter part of Appendix C, pages 45-46 of this report. As for 
the sample of size n = 100, see the passage from [9] quoted in connection with (C7) in 
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this report. The point involved here is that, although the conditional probabilities p, 
in the WB method cannot be assumed to be normally distributed, the means p of 
samples of size n are approximately normally distributed, by the central limit theorem, 
q and the larger the value of 1, the better is the approximation to normality. Reference ag 
(9, p. 1111 states that this normality assumption is quite accurate in most cases 
for n > 10. Hence, we can assume that, in this application, it is very accurate for 
n= 100. and that the failure of the 50% confidence interval to include the apparentiv ’ 
correct value P, = 0.044321 cannot be attributed to the fact that the means § in the 
WB method are not from an exactly normal distribution. It should be borne in mind 
that, no matter how large the sample size n may be, the probability that the 50% 
i confidence interval wiil cover the true population mean, P, or Po» ds, by definition, , 
a only 0.50, For higher confidence levels such as 90%, the confidence intervals are 
4 longer than for the 50% level. 


a ! This last-mentioned calculation (WB method with n = 100) used something of 
f the order of 36 seconds on the CDC 6700 computer, since this and a similar run, 

using 100 impact points on a different case, took a total of 73 seconds. If fewer 
simulated impact points are used, the machine time is proportionately shorter. But q i 
rt: since the P.E, (probable error) is approximately proportional to |/n, where n is the 
4 number of simulated bomb impact points, the smaller valucs of n are associated with 
“i relatively large P.E.'s, 


ht. The computing time per case by the DJH method, using double integration , 
§ (2-12) or (2-13) and subroutine PKILL, was between 1.0 and 1.5 seconds for those : 
a cases which were timed, of the 22 cases referred to above, in comparing the DJH and | 
g j WB methods. However, the computing time per case was much larger, in some cases a 
q | of the order of one minute or more, in some of the extreme cases which were com- 
puted in selecting a sct of input values for the tables in this technical report. This led ; 
{ : to some changes in the double integration subroutine used, in an effort to increase rt 
| 


computing speed while maintaining acceptable accuracy, This point is discussed in 
some detail below. 


Simple tests were included in the program for detecting some, but not all, of the 
cases in which the result will be P, < 0.005 and, therefore, P, will be given as 0.00 g 
in a two-place table, and for detecting some of the cases resulting in P. > 0.995 (and, . A 
therefore, P, = 1.00 correct to two decimal places). We will give a discussion in terms OE 
of a small but unspecified e, although the actual tolerance used in computing Table 4 
E3 was e = 0.005, so that this table may be regarded us correct to two decimal | i) 
digits. If the values in a table ure subject to crrors not exceeding 0.005, the table 
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is conventionally regarded as coirect to two decimal digits, although there are excep- 
tions, as for instance when the computed value is 0,824, rounding off to 0.82, with an 
error of 0.003, making the true value 0.827, which rounds off to 0.83 rather than 0.82, 
More generally, if the errors do not exceed § x 107°"! or (1/2) x 107", the table 
is conventionally regarded as correct to n decimal digits, although it would be more 
correct to state that the error, after rounding, does not exceed one unit in the nth 
decimal place. 

For the coordinates of the points (his Kaun ) involved in 
these tests, and further details, see Appendix D. 


) and (h 


max? K ax 


In the first of these tests, for the detection of cases in which P, < ¢ and P, 
will be given as zero to the appropriate number of digits, a point (h,, in? Kinin) with 
non-negative coordinates is determined, such that, for every point (h, k) on or inside 
the ellipse, we have |h| > Hain’ \k| > Kin’ Consequently P(h, k,o,,0,>R)$ PCR in? 
Kinin! °n? 2%» R), since the elliptical c coverage function decreases with increasing jhi and 
also with increasing |k|, Further, since P, us given by (2-12) or (2-13) is simply an 
average value of P(h, k, ¢,, 0,, R) over the interior of the ellipse, by (2-17), it will be 
true in such a case that Py < PChi is Kins Ono 2,» R) with strict inequality, since we 
always have strict inpquallty at interior points of the ellipse. Hence, if P(h, i. Kia» 
0,9» R) < «, requiring only one call to the PKILL subroutine for the selected values 
of h and k, we must have P, < e, This test will fail in some cases in its objective of 
saving computing time, since it ls possible to have Py <e< P(tains K eis “1 a,» R); 
but it cannot lead to an incorrect decision that P, < e, since, if P(h 
o,, R) § e, we must have P, < «. 


min’ Kinin? Oy? 


We remark here that, in the preliminary computation of some of the pages of the 
present main table with the coarse tolerance of e = 0.050, the criterion Phin? K ae 
0,29.) R) $ € was replaced by P(t K min? 2nt Px? R) § €/2, but this is “on the safe 
side” and certainly not less rigorous than the standard criterion involving e. The objec- 
tive was to get a smoother decrease in the values of P, with increasing u and v, as 
illustrated by the following example. With « = 0.050, suppose that P( hl. Kuins Op) 
a,» R) = 0.045 and P, = 0.040 if computed by the regular DJH program. If the 
regular e criterion is applied, P, = 0.0 (with accuracy to only one decimal digit if 
e = 0.050). But PCH,» k min? 7h? Pk? R) and the computed P,, are both greater than 
e/2 ™ 0,025, so that P, is ‘isk printed as 0.0 under the e/2 criterion, The value 
P, * 0.040, though seobably not accurate to three digits, may have a substantially 
sinaller actual error than the maximum error indicated by the analysis, and the value 
0.040 contributes more than the alternative value 0.0 to a smooth decrease in P, 
with increasing u and v. 
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In the analogous test for E. > 1 — e, a point hs 55 K vax? is found, such that, 
for every point (h,k) on or inside the ellipse, we have {hj < Nex? kl <¢ Kinax’ Then 
PR, > P(h., wa Regs 0,5 0%, R) (by reasoning similar to that in the corresponding step 
in discussing the criterion for P, <e). If a call to the PKILL subroutine determines 
that PCH» Knax? hs %» R)> 1 — ¢, we must have P, > 1 —e. Here it is possible 
to have Ph > 1 —¢> Pthi kK 9,1 9%» R), in which case the test fails'to save 
computing time; but if PC. ax? Kmax? [hr Ox» R) > 1 —e, we must have Bol =e 
and, therefore, the result to be printed in the main table will be P, = 1.00 «(to 
the appropriate number of decimal digits) without the necessity of evaluating the 
double integral by the regular DJH program. 


Originally, it was planned to compute a table of 16,384 values of P,, with R 
normalized at 1 and the remaining input parameters as follows: 


a= 1,2.8,10 o,/o, = 0.1,0.5,0.8,1 
b/a = 0.1,0,2,0.5,1 u = 0,2,5,10 

@ = -45°,0°,45°,90° v = 0,2,4,6 

0, * 0.5,1,2,5 R=] 


If all combinations are taken, the resulting table has 4” or 16,384 values. 


However, it was found from theoretical considerations involving somewhat more 
than the test described above for Py > 0.9995, that such a table would contain no 
values exceeding this number (an unrealistically fine tolerance of e = 0.0005 was under 
consideration at this time). A key fact here is that the ratio of R (=1) to the smallest 
value of a, (0.5) is only 2. 


Several skeleton tables, cach of 128 values, were computed with the object of 
selecting a set of input parameters which would result in a better balanced table in the 
sense of having a greater proportion of values close to 1.000, and also a reasonable 
proportion close to 0.000. The following is a typical set resulting in 128 cases, with 
R = | in all cases; 


0, = 0.1,2 u = 0,10 
o,/o, = 0.1,1 v = 0,4 
a= 0.1,8 @ = 459,909 
b/a = 0.1,1 R= 1 
18 
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The short tests described above for detecting cases in which R< 0.005 or P,. > 0.995 
were used, with positive results in a substantial proportion of cases in these skeleton 
tables, saving computer time. 


Later it was decided to normalize by taking o, = 1 (rather than R = 1) with 
5, $ 0,.always, A number of sets of 128 test cases were computed with this normal- 
ization, a typical set being: 


a, 71 u = 2,0,3.0 
o, = 0,2,0.5 v= 1,0,1.5 

a= 0.5,2.0 “@ = 0°,45° 
b/a = 0.5,1.0 R = 0.5,1.0 


However, in a number of these test cases, the computing time on the CDC 6700 
machine was unexpectedly long, sometimes more than one minute per case. It was 
concluded that the numerical double integration subroutine was probably using an 
unnecessarily fine mesh in many of these cases, and study was given to methods of 
avoiding this. The PKILL subroutine is a relatively slow ono, requiring something of 
the order of 5 milliseconds even for 3-place accuracy on the 6700 machine, The 
CIRCV subroutine (circular coverage function) is an order of magnitude faster, but 
can be used if, and only if, o, = o,, as has been stated. At 5 milliseconds per case, 
only about 200 elliptical coverage functions per second can be computed to three- 
place accuracy. Hence, if the double integration subroutine should use a mesh of 128 
subintervals or 129 points on each dimension of the rectangle of integration, there 
would be 129? or 16,641 calis to the PKILL subroutine, requiring more than 83 
seconds in calls to this subroutine at 200 calls per second. 


The error analysis, resulting in improved speed in computing the double integral 
of (2-12) or (2-13), while maintaining acceptable accuracy, is given later in this section. 


A number of test runs of 128 cases at a time, with o,, always having the value 1, 
were made. Also, curves were computed and plotted in the hk-plane for P(h, k, o,, 
o,, R) = 0.005, and 0.995 foro, = 1, and various values of co, and R. This facilitates 
the determination of ranges of values of the input variables in the double integral of 
(2-12) which would result in values of P, = 0.000 or 1.000 to three decimal digits, 
by the tests described earlier in this section. On the basis of all these tests, it was 
tentatively decided to compute a table of 16 pages with 1,024 values to a page, mak- 
ing a total of 16,384 (= 2'*) cases, with the following input values: 


ay 
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g, * 0.2,0.5,0.8,1.0 (4 values) 
u = 0,0.5,1.0,3.0 (4 values) 
v = 0,0,25,0.50,1.50 (4 values) 
¢@ = 09,459,90°,-45° = (4 values) 


a = 0,1,0.5,1,0,2.0 (4 values) 
b/a = 0.2,0.5,0.8,1.0 (4 values) 
R = 0,5,1,0,2.0,4.0 - (4 values) 


Taking all combinations, the total number of cases is 47 = 2'4 = 16,384 as stated 
above, 


Later, it was decided to include tho value R = 0.25 on the ground that this 
might be a more realistic value of the lethal radius R than the larger valucs, from the 
practical military point of view. This would make a total of 20 pages of tables, or 
5 x 4° = 20,480 cases. 


When the computing of tho’pages was started, it was found that some of the 
pages had an excessivo number of values near 0 or 1, and additional pages were com- 
puted by modifying and augmenting the set of values of the Input variables, especially 
u, on some of the pages, The resulting total for Table E3 was 51 pages, containing 
59,904 values uf P,., as stated earlier in this section, page 11. 


A tolerance « = 0.005 was used on all pages except 61 and 65, on which the 
tolerance wis 0.001 (see page 12 for further comments on this). The computing 
time per case varied widely, being 0.10 second or less on some pages, and 0.85 
second or more for other pages, The overall average, from available recorded times, 
was about 0.4 second per case. 


In using the DJH program, based on (2-12), arbitrary positive values can be given 
to o,,9,, 8, b (4> b), and R, and arbitrary real values to u, v, and --909 < g< 90°, 
There is nothing in the analysis which rules out the case o, > 9, Ifuor vor both are 
negative, the corresponding positive value or values could be taken, with appropriate 
changes in ¢;for example, u = —4, v = 2,¢ = 30° is equivalent tou =4, \ = 2,¢ = ~30°, 
as is clear if a figure is drawn, the elliptical coverage function being syrametric with 
respect to both axes. 
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But if the table in this report is used, with or without interpolation, it must 
be noted that, in the table, 0, = I and o, > 0,. Hence, it may be necessary to nor- 
malize the input variables (except @) by making o, = 1, and ifo, > 9, originally, the 
roles of h and k must be interchanged, If we have the case o, =5,0, =2,u=3,v= 2; 
a= 2, b = 1.5, @ = 30°, R = 3, we simply divide all distances by 5 and obtain as an 
equivalent case 0, = 1,0, = 0.4, u = 0.6,v = 0.4, a= 0.4, b =0.3, ¢ = 30° (no change 
in 6), R = 0.6. 


Now, suppose we have a case with o, > o,, say o, = 2,0, =5,u=3,v=2, 
a=2,b=1.5,9 =30°,R=3. We first divide all distances by the larger o, 0, in this 
case, and get the cquivalent case o, = 04,0, =1,u=06,v=04, a= 0.4, b= 0.3, 
» = 30°, R = 0.6. Then, reversing the roles of h and k, we will also have to reverse 
those of u and v and make an appropriate change in @, obtaining the equivalent case 
o, 1,0, =04,u5 0.4,v=0.6,a5 0.4, b = 0.3, ¢ = 60°, R =0.6. This should be 
clear from a study of Figure 2, bearing in mind the symmetry of the elliptical coverage 
function with respect to both axes, Note that ¢ = 60° and not —30°, 
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Figure 2. Interchange of o,, 0, in Example 


Here we discuss the numerical evaluation of the double integral (Figure 3) 
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by a method which has been chosen as an efficient one, using two tolerances ¢, and 


e, and an overall tolerance ¢. 
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Figure 3. Area of Integration for Double Integral 


' These three e’s may be given appropriate positive values depending on the accuracy | " 
required; for example, the overall tolerance e, on which ey and ¢, depend in the pre- | 
sent analysis (see (3-15)), is sometimes given a value of 0.01 in crude calculations, and 
sometimes a value of 0.00001 if much more accurate results are desired. | 

| 


In (2-12) and (2-13), the variables of integration are r and 0. We leave open the | 2 
possibility of identifying r with x and 6 with y, or the reverse. ; 4 


We regard the integral as an iterated one, 


E d b d 
r= f | fi Fo oso fin ey (3-2) : 
C a c ' 


b 
For a fixed y, the integral f(y) = I F(x, y) dx is approximated by Simpson’s rule, 


bui the number of subintervals for two different ordinates or values of y is not 
necessarily the same, Hence, we do not have a two-dimensional grid of the less | 
efficient type sometimes considered in the literature, for example in [7, p. ‘186, | 
problem 54}, 


ts For y = for arbitrary s on (c, d), let 


- é b 
re , 7 exact value of f(y) = f(s) =f F(x, s) dx (3-3) 
a 
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1, = Simpson’s rule approximation for I, (3-4) 


in general, a bar will indicate an exact value, as in E, I in this case indicates an 
approximation. 


1, #1, +5, (3-5) 


where E, is the error committed by using the Simpson’s rule approximation with the 
chosen number of subintervals of [a, b]. 


For a given s, ie., a given value of y, let a sufficiently large number of subinter- 
vals of [a, b] be chosen so that 


I, - 1,1 = IEl<e, (3-6) 
where e, is one of our fundamental tolerances. 


We are really comparing successive approximations, say (I,), and (I,)», 
with each other, where n or 2n denotes the number of subintervals into which 
[a, b] is divided, rather than comparing any approximation with the exact I, 
since in practice we do not know the value of the latter. We do not in this 
analysis derive rigorous error bounds for integration by Simpson’s rule. 
Numerical analysts have developed several such rigorous bounds, using the 
fourth derivative of the integrand [7, p. 75, (3.5.23)], or evaluating a line 
integral in the complex plane, or using other methods.’ The calculation of any 
of these rigorous bounds would require a relatively large amount of machine 
time, and, moreover, such bounds are often ultra-conservative, i.e, much larger 
than the actual error committed. It is felt by the present authors that 
sufficient accuracy for practical computations can be attained by developing 
simple tests involving convergence of successive Simpson approximations to a 
common limit. Such tests can give a high degree of confidence in the reliability 
of the results, thoigh not absolute certainty. We may decide that (3-6) is 
satisfied if two successive iterates agree with each other within ¢,; or, if we 
are more conservative, within €y/2; or, of three successive iterates, say (I,),, 
(I,)o, and (1,),,, the largest and smallest are within €, of each other. We 
chose the first of these. We state in (3-17) below, how the tolerances are 
handled in the present NSWC program, but the program is flexible in that 
changes in the tolerances are easily made. 
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Now by definition, 
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= exact value of double integral (3-7) 


Also, I, is a well-defined function of s or y once a numerical value of ¢, has been 


assigned and a criterion for the satisfaction of (3-6) has been chosen as discussed 
above. 


d d 
We assume also that I I, dy or i I dy exists. I, or fore, y) dx (3-3), the 


exact value of the integral from a to b, is a continuous function of y, by the conti- 
nuity of the integrand of (2-12) as a function of two variables. A consideration of 
the manner of forming the Simpson sums shows that the approximation I, is also a 
continuous function of y in any interval on the y-axis in which the number of sub- 


divisions of [a. b] in forming the Simpson approximation to Sor, y) dx dovs tiot 
change. But l, does have a discontinuity, in general, at a Boe where we switch 
from 16 hy 32 aunienaie for example, in approximating hg F(x, y) dx. Suppose 

=> F(x, ¥,) dx uses 16 subintervals for some y = y,, but J. F(x, y) dx uses 


2 subintervals for some neighboring values of y, say for all y in the open interval 
(Yo: ¥o + k) for some positive k. Then I, as a function of y, will in general have a 
discontinuity for y = y,. Nevertholess, the integrand in (2-12) is such a smooth 
slowly changing function of its arguments that we can be confident on heuristic 


d : va 
grounds that the integrand in t I dy has a small finite number of discontinuities. 


“d 
But the Riemann integral I I dy exists if, and only if, the set of discontinuities is 
of Lesbesgue measure zero [8, p. 230]. This includes the present case, in which the 
-d 
discontinuities are assumed to be finite in number. Hence, the integral tI ly dy or 
d 
I I, dy exists. It is shown below in (3-8) that it differs by less than (d -- c) €, 
( we 
from the exact value of the double integral, J , dy orl, ie, it approximates T 
= n 
arbitrarily closely for small ¢,. Also, the Simpson sums de Zo k, I, which are 
used in approximating the integral with respect to y, the k,’s being the Sitapson 
d 
multipliers 1,4,2,...,4,1, converge ‘o the value of the integral te ly dy as the number 


24 


oa aly ge ys epbsana dee ov a see 
jini ing ; 


: of subdivisions of [c, d] becomes infinite (see (3-10)). It is not difficult to prove that 
| the Simpson sums, assuming that they converge to the true value of the integral in the 
; integration of a continuous function over a finite interval, also converge to the value 
t i? of the integral if the function has a finite number of discontinuities with a finite jump 
\ at each such point. 


i 

Assuming then that Ps I, dy exists, and assuming that (3-6) holds everywhere on 
} {c, d], ie., that we take a sufficient number of subintervals of [a, b) so that (3-6) 
holds for all y on [c, d], we have 


| [fiw fal [foveal L 
s fi-4| av< fie aye - oe, (3-8) 


or 


: r- fi, ay| <@-oe, (39) | 
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Now, suppose that we take a sufficiently large number of subintervals on [0, d] 
so that 


d—c 
vo dy ae k, 1. <e, (3-10) 


where the sun uses the | 3 for the finite number of points taken on [c, d] (=n + 1 
points where n is the number of subintervals of [c,d])., Here, e, is another funda- 
mental tolerance, like e,. As in the case of e,, we can require two successive iterates 


is zk, 7 to agree with each other within e,, or within e, /2, or choose some other 
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some arbitrary criterion for the satisfaction of (3-10). 
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4 criterion. We do not know the value of the integral J I, dy, so we are forced to adopt | 
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With (3-6) and (3-10) satisfied, we now have 
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correct within e. From (3-13) and (3-14), we equate e with (d — c)e, +e,, so that 
we have 


jr -S=* sy, <(d—c)e, te, =e (3-15) 


~- d-c d—c 

Ir- a zk I, = [t- fi,ays fi, eer k, | (3-11) 
d—c d—c % 
| [r-$== si, 1, | $ |r fi, as| + [fiey-S2 | (3-12) 
, 1 
: - d-c I j 
ar ree 4, <(d —c)e, te, (3-13) ; 
Now suppose that ¢ is some overall tolerance, and we wish to choose ¢, ande, | 
; in such a way that | 
[r- 2 x 1, | <¢; (3-14) 
3n i 
: that is, we wish the computed Simpson sum approximating the ¢uuble integral to be | 

| 


Equation (3-15) clearly allows some latitude in deciding on the individual values 
of e, and ¢,. The present NSWC/DL program takes ty. = eis from which follows 


[l+(d ~ 0] «, (3-16) 


TO me me eee 


= € 


¢, Ze, * T*a@oo (3-17) 
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If we wish to experiment with various values of the ratio k =e, /e,, a simple generaliza- 
tion would be, from (3-15), 


(d-c)e, tke, we (3-18) 


ere es! 
° k+(d —-c) 


| 
(3-19a) 
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1 (3-19b) 
where k is any positive real number, less than, equal to, or greater than 1. The case 
k = | reduces to (3-17). It is possible that computational experience would show 
that some value such as k = 2 or 1/2 is generally more efficient than k = 1. The 
determination of such an optimal value, or approximately optimal value, of k, would 
be purely an empirical matter. We could try various values of k fora fixed value of e, 
and regard as optimal that valuc of k which resulted in the smallest expenditure of 
machine time while maintaining the specified accuracy. 


A simple modification of this error analysis is necessary in case the constant 
factor 1/(ra?) outside the integral signs in (2-12) and (2-13) is not included as part 
of the integrand at every step of the analysis, but only once at the end in computing 
the true value of P,. This is the cuse in the DJH program. 


Let the integrand in (2-12) and (2-13), not including the factor 1/(7a?), be 
denoted by F*(0, r), and, in general, let the asterisk or “star” denote this omission 
of the constant factor, Then (2-13) can obviously be written in either of the forms 


2 P I * ("0,1 a0] d 
Gi = wT. r 
K Jo |¥o 


ON I A ie es ES OO reli me ee, on a 


i } 1 a] p2n 

k : el, f F*(o, 1) de | dr (3-20) 
4 t ra? JO {JO 1 4 
A i a 

y i where ! 
Hf t 4 
' F*(9, 1) * P(h, k,0,, 0,, Rt (3-21) ’ 
F(, r) = fh F*(6, r) (3-22) q 
na? x 


and h, k are given by (2-11). 
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Also, 


wm ff F* do dr (3-23) 
and 
P =[ : I* 3-24 : 
K aa (3- ) 


Now, following the foregoing analysis but using the “star” in the sense which 
: has been described, we evaluate the double integral I* to an accuracy such that 


jt 228 


Bk, Ik] <et (3-25) 


(with e* to be determined in terms of the overall tolerance ¢), with 


eet 3-26 
‘0 1 Td ake 


F as in (3-17). We will then have 


! 
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3 a (Te Sof yk ye] < Set 3-27 
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y 
Lf 
‘ But we wish to have 
i 
| 
; and so we put e*/(ra*) = e or ! 
e* = ma? + ¢ (3-29) | 


st a 


tase om 


and therefore, from (3-26), 


ef = en = rls cea Ree 
1+ (d—-¢) 1+(d —c) 


The tolerances are handled in this way in the DJH program. The constant ¢ is set at 
0.005 in the present program, but can easily be changed to some other desired value 
such as 0.001 or 0.05. 


We close this section with a note regarding algebraic signs of h and k and limita- 
tions of subroutine PKILL. 


A glance at the integrand of P(h, k, o,.0,, R) (see page 3) shows that P(h, k, on) 
a,» R) is symmetric with respect to both h and k. Consequently, the subroutine 
PKILL takes input values of h and k, which may be negative, and replaces them with 
{ht and {kl respectively. Moreover, the PKILL subroutine has the limitations that 
Ihl/o, © 600, IkI /o, < 600, Although it is unlikely that such extreme values will ever 
occur in the integration of P,, the precaution is taken to sense on these inequalities 
with an immediate exit from the routine if they are not satisfied. By means of 
Equations (D | b de g), page 51, a rectangle containing the entire ellipse of Figure 1, 
page 3, is determined. If this rectangle is contained within a rectangle, centered at the 
origin, of width 1200 o, and height 1200 o,, the inequalities Ihi/o,, © 600, IkI/o, < 
600 must be satisfied everywhere in the field of integration, If not, a value P,, = —1 is 
returned by the program, 


The standard deviations 0, and 0, are subject to the condition 
Wis < g,,/o, «< 15, (2, page 26]. 


If the condition is not satisfied, the program exits and also returns a value of 
P=~1, 
k 


| 
i 
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SPECIAL CASES 


Two special cases are worthy of mention here: (a) a, 79,3 (b) a=b, ie, 
the ellipse is a circle. 


(a) We have assumed that o, = o,, with strict inequality in general. But 
ifo, * 9, * 4, the distribution is circular normal, and if we divide all coordinates by 
a, we have a standard circular normal distribution, ie., with o = 1. We would write 
h = h/o, k= k/o or replace h by oh, and similarly replace k by ok. Also, we would 
divide, u, v, a, b, and R by o. The angle would be unchanged by the transformation. 
Ifc, =o, =o = 1, the case is simplified and the original symbols, h, k, u, v, a, b, 
and R can be retained. 


We would then evaluate a double integral as in (2-12) or (2-13), but call the sub- 
routine CIRCV, computing the circular coverage function P(R, d), where d = +/h? + k?, 
instead of subroutine PKILL, which computes the elliptical coverage function Ph, k, 
Ons Op, R). The former subroutine is at least ten times as fast as the latter subroutine. 


This substitution of subroutine CIRCV for PKILL, and the reduction dis- 
cussed below of the double integral to a single integral, also using CIRCV, can be 
done if and only if o, =0,. The reason is that if o, #0, a transformation must be 
made, in order to produce a circular normal distribution, in which unequal factors 
of stretching or compression along the h- and k-axes are used. This would convert 
the lethal circle of radius R into a lethal ellipse and this is contrary to the assumption 
that the bomb’s lethal area is a circle. 


If o, = 9, = 0, involving an initial similarity transformation if o * 1, as dis- 
cussed above, we can even express P, in ternis of a single rather than a double 
integral, but at the cost of solving a fourth degree equation in general at every inte- 
gration step to find the points of intersection, possibly four in number, of the circle 
h? + k? =r? (whero r is the variable of integration) with the ellipse. 


Figure Al shows a case in which a circle with center at C and radius r (also a 
circle of radius r + dr) intersects the ellipse in four points. The arcs of the first circle 
(radius r) subtend central angles @, and 6,. Let 9(r) = 6, +6,. Then the sum of 
the two shaded areas is (0, + 4,) r dr or 6(r) r dr neglecting infinitesimals of higher 
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Figure Al, Case o, = 0,, Ellipse Intersected at Four Points 


order than the first. This area can be described as the area occupied by points that 


are in the interior of the ellipse and at a distance between r and r + dr from the center 
C of the distribution. Because of the uniform distribution in the ellipse of possible 
positions of the point target, the probability that the target lies in the shaded area is 
6(r) x dr/(rab), i.e., the ratio of the shaded area to the total area of the ellipse. All of 
these points in the shaded area are at distance r from C, neglecting infinitesimals of 
higher order. Hence, given that the target lies somewhere in the shaded area, the 
conditional probability of its destruction is P(R, r), the circular coverage function as 
computed by subroutine CIRCV. Hence, the total probability that the target lies in 
the shaded area (for a given r) and is destroyed is 9(r° r dr P(R, r)/(rab), and the overall 
kill probability P| is computed by integrating this expression from r, tor, , which are 
respectively the minimum and maximum of values of r such that a circle of radius r, 
center at C, contains points that are on or inside the ellipse. If the point C is on 
or inside the cllipse, r, > r, = 0, but if C is exterior to the ellipse, r, > r, > 0. Thus, 
we evaluate the single integral: 


f 
Py = [0 P(R, t) dr (Al) 
nabJr, 
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The determination of r,, r, and 6(r) in general involves the numerical solution of 
fourth degree algebraic equations, but is simpler in special cases as discussed below. 
If C is in the interior of the ellipse, so that r, = 0, we take 6(r) = 2n for all values of r 
such that the circle of radius r (center at C) lies entirely in the interior of the ellipse. 
If C is on the ellipse, we also have r, = 0, but for all circles of radius r > 0, 6(r) < 2n. 


A special case which was programmed and used in checking out the general pro- 
gram of this report is shown in Figure A2, In this case, the limits r, and r, of (Ai) 
are obvious by inspection, and the computation of the function ar) of (AL) involves 
only square roots and not the solution of more general fourth degree equations. 


MONA G, 
Cp , aN, 


Figure A2. Special Case Programmed 
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Here, we have (> oS 1, u = 0, 6 = 90°, v > a> b, so that the point C is 
exterior to the ellipse. From the figure we clearly have ry =v-—a,r, =vta. It is 
intuitively clear that for v - a <r < v +a, the circle of radius r, center at C, will 
intersect the ellipse in two, and only two points, with equal ordinates by symmetry, 
and this is confirmed by analysis. The common ordinate y of these two points of 
intersection is given by: 


Va? — 67) ~~ h2) 4 bev? __ h2 
a (a oe b?v (A2) 
and, when y has been computed, the abscissas of the two points of intersection (circle 
and ellipse) ste + ys? — y?. Hence, the angle @ (Figure A2) is arc tan (¥ 1? — y?/y), 
and the entire central angle corresponding to the shaded area is twice this or 2 arc 
tan (vr? — y?/y). Therefore, using (Al), the kill probability P, is given by the 
integral: 


i arc tan P(R, r)r dr (A3) 


nab JV--a 


Ps 2 vta fi = y? 
y 


with y at each integration step, or each value of r, given by (A2). 


Although a program based on Figure A2 and (A2) and (A3), with a> b, was used 
in checking out the general method of this report, we would have a still simpler spscial 
case if the ellipse were a circle of radius a, with a < v (see Figure A2 and imagine that 
a= b). In this case, the ordinate y of the points of intersection of the two circles, 
h? +k? =? and. h? + (k - v)? = a?, is given by: 


r2 + y? -- a? 
y a ae (A4) 


and P,, would be given by (A3), with y computed by (A4) instead of (A2). 


(b) If u = b, implying that c = 1 and the ellipse is a circle (in the general case 
where ¢ and o, are not necessarily equal), the positions of the axes of the ellipse 
are indetcrminate, and for convenience we can take ¢ = 0 so that the x- and y-axes 
ure parallel to the h- and k-axes respectively. This would simplify (2-3) which would 
become 
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Also, c would be | in (2-5). If o, # o,, we would still have to compute a double ae j 
integral as in (2-12), using subroutine PKILL, but with the simplifications which have 
been noted in the coordinate transformations, (2-3) and (2-11). 
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ALTERNATIVE METHOD OF COMPUTING P,. 


We outline here an alternative method of computing the kill probability P,, 
this thme using true polar coordinates t, @ in the ellipse. The curvilinear coordinates 
introduced in (2-5) and also denoted as r, 9 are not polar coordinates in the ellipse 
ifa> b. 


Figure B1 shows a quadrani of the ellipse (see also Figure 1 in Section 2), 
First, we wish to find the coordinate r for a point such as E on the ellipse, with 
given polar angle @, to establish limits on 1 in terms of 9. The equations of the line 
OE, and of the ellipse are: 


sweet pee ne an 


y™ x tan 0 (B1) 
2 2 
. foe (B2) 


i Figure B1. Polar Coordinates in the Ellipse 


i Solving these equations simultaneously, we find that the rectangular coordinates 
{ (x, y) of the point E are: 


x = (ab/S) cns 6, y = (ab/S) sin @ (B3) 


} 
i 
y 


and, therefore 


r= Vx? + y? = ab/S (B4) 


where 


S = Va? sin?0 + b? costa (BS) 
Equation (B4) provides the upper limit on the variable r in the subsequent double 
integration. 


Now, considering a polar element of area r dr dé (Figure BJ) at a point (r, 6) 
inside the ellipse, we must, as before convert its rectangular coordinates (r cos 6, 
r sin @) to the Chk system, by (2-3). The result is: 


h*u+rcos@ cos¢ —rsin@ sing #utrcos(d +¢) 


a 
k™v+rcos@ sing +rsind cosg =vtrsin(6 +o) (B6) 


By reasoning similar to that used in setting up (2-12), we obtain the following double 
integral in terms of the present true polar coordinates r, 9 


P= 1 2m ake P 9+ +rsin(o + R] r dr dé (B7)* 
y mae A f [utr cos(a +), vtrs $),0,,0,,R) rdr 


The P functions, as in the case of (2-12), are computed by subroutine PKILL. 


“See note, page 29, regarding h and k. 
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PROBABLE ERROR IN THE WEIDMAN-BRUNNER METHOD 


In the WB (Weidman-Brunner) method as described in [1] and in the present 
technical report, the standard statistical method for determining the probable error 
and confidence intervals is one using the Student t distribution, [4, p. 224], [6, 
p, 178], and [10, p. 189]. For a given problem, o,, o,, R, u, v, a, b, and  speci- 
fied, there exists a true kill probability, P, or p,, which we are estimating or approxi- 
mating by the sample mean, f, of a sample of size n (the number of simulated bombs), 
For a fixed and sufficiently large n, these sample means { B} are from an approximately 
normal population with mean P, or p, and standard deviation o (for individual prob- 
abilities p,), the standard deviation of the sample mean § being of /n. The point 
which makes the t distribution the appropriate one to use in computing confidence 
intervals is that o (as well as P, or p,) is unknown. The method cited from [4], 
(6), [10], and given in many other statistical texts, is designed for this situation, 
finding confidence intervals for the mean of a normal population with unknown 
standard deviation o. 


However, in the present technical report, we do not develop the above-mentioned 
method in detail, since it is not used in {1] or in the Weidman-Brunner computer 
program. Instead, we follow this computer program and [1] more closely, and, with 
each sample, estimate the population standard deviation o by means of the sample 
standard deviation. It is believed that this method if followed up correctly will give 
results that are roughly equivalent to those given by the method which uses the t 
distritution. 


The method used in [1] for computing the probable error (P.E.) is open to 
criticism in several respects. A study of the FORTRAN program given in this reference 
shows that the P.E. is given by the formula: 


ft 1/2 
P.B. = 0.6745 | S° (p, ~ 5)? | (C1) 


imi 


where n is the number of simulated bombs, p, is the i-th conditional probability 
determined as cxplained in Section 1, and B is the mean of the n p,’s,ic., B isthe — 
approximation to the kill probability P,. The integer n is usually taken as 25 in the 


WB method, although the values n = 50 and n = 100 were used in some of the cases dis- 
cussed in Section 3 of the present report. The radical in (C1) represents the standard 
deviation o of the sample of size n, and the constant 0.6745 is a familiar one in 
statistics, where P.E. = 0.67450, since, for a normal distribution with mean O and 
standard cieviation o, we have: 


0.67450 j x? 
2 


1 
eens ex - 
of 2n Joe , 


J] dx = 1/2 (C2) 
20 
or, half of the items lie between the limits + 0.67450, 

The denominator under the radical in (C1) should be n — 1 instead of n in order 
to get an unbiased estimate of the standard deviation [6, p. 154); actually, the use of 
n ~ 1 gives an unbiased estimate of the variance a?, and a very small bias in the value 
of o itself. For large n, there is no significant difference between using n and n ~ 1 
in (C1). 


A second and inore serious criticism is that even if we substitute n — 1 for n in 
(C1), the result is the P.E, for an individual item p, and not for the sample mean p. 
To compute the P.E, for the latter, an extra factor of 1 is needed in the denomina- 
tor of (C1) as shown by the following discussion. 


The conditional probabilities p,, p,, .., p,, in the WB method are independent 
random variables with identical (though in general unknown) distributions. The dis- 
tributions are identical since, for each of the n simulated bomb drops, the procedure 
is oxactly the same. But we cannot assume that this common distribution Is normal, 


since it involves a number of factors including the uniform distribution of the point 


target in the ellipse. 


In this situation, we can get a good approximation for large n, for the variance 
(hence, for the probable error P.E.) of the sample mean p by means of the central 
limit theorem [4, p. 136] and [6, p. 191). Reference [6] states (changing the nota- 
tion slightly to agree with our present notation) that if p,,p,,..., p, are independent, 
identically distributed random variables, each with mean p, and variance o?, and if 
y, =P, +P, ++ P, ™ np, and 


Y, ~ ™Py_ n(p—p,) n(B~p,) . 
ey a (C3) 


ovyn 


| then Z, is “asymptotically normal” with mean O and variance 1, in the sense that the 
distribution function F,(p) of Z, satisfies the relation 


ase 


1 P 
lim F Smee —u? /2 y 


Here, p, is the true kill probability for a given problem of which 6 is an estimate on 
the basis of one sample of n bombs. We think of a large number of samples of size n. 
The random variable p, takes a large number of values in these successive samples, 
with mean p, and variance o? and similarly with Pas Py, ete, | 


; Since Vn (B — p,)/o has approximately mean O and variance 1, as stated in the i 

formulation of the central limit theorem in [6, p. 191], from which (C3) and (C4) are : 

; taken, 6 — p, must have mean 0 and variance o?/n and, therefore, f, the sample 

{ mean, must have, asymptotically, mean p, (the true kill probability) and variance 
o7/n, establishing the fact that the variance of the mean is 1/n times the variance of ; 

\ an individual p,. 


The facts that the mean of all these sample means is p,, and the variance is | 
a4/n, could have been proved by simpler theorems, [6, pp. 135 and 137] but to | 
prove that the sample means p are asymptotically normally distributed (and this fact 
is needed in computing the probable error) requires the central limit theorem. Also, | 

; if the individual random variables p, were known to be normally distributed, we could 

i get very simple and exact, not merely asymptotic, results for the distribution of the 

: sample means p. But under the present circumstances, with the P, ‘s not normally i 
| distributed, we apparently cannot do better than the present asymptotic analysis, 
y 


Since the variance of the sample means § is asymptotically o?/n, as has been 
shown, the corresponding standard deviation is o//n, where o is the standard 
deviation for an individual value Py For this o we have to use the estimate { 
f [s(p, - BP n ~ 1))'/2, square root of an unbiased estimate of the variance as 
% stated above. We multiply this by Jn ta get an estimate for the standard devia- 
i tion of the sample mean f, and by 0.6745 to get the probable error, P.E. (see (C2)). 

Hence, we have 
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P.E. = 0.6745 [(p, — B)?/n (n ~ 19] Mig (C5) 


Let the P.E. of [1] be denoted by P.E.,, (see (C1)). Comparing (C1) and (C5), we 
see that P.E./P.E.,, * 1/y n-~-1,or 


PE, = ——Yh. (C6) 


P.E.,,, i8 the value printed out for the probable ertor if the FORTRAN program of 
[1] is run on a computer. These printed values were corrected as indicated in (C6) 
in interpreting the results of computed cases. 


Although the distribution of the sample mean f in the WB method approaches 
normality as the sample size n becomes infinite, there appears to be no simple means 
of determining how fast its approach to normality is, or how large the sample size 
should be to obtuin a meaningful estimate of the probable error in the computed 
approximation p to the true kill probability p,. In (9, p. LIL], in a-discussion of 
this application of the central limit theorem, the following statements are made: 
“In practical terms, a normality assumption for the sampling distribution of X 
becomes reasonable in many cases for N > 4 and quite accurate in most cases for 
N > 10. Hence, for reasonably large sample sizes (4.34) applies to the sampling 
distribution of X computed for any random variable x, regardless of its probability 
distribution function.” Equation (4.34) [9] is 


Ku) VNL, (C7) 
qa 


* 


which is stated to have a standardized normal distribution if x is normally distributed, 
The normality is only asymptotic if x has some other distribution. The expression 
on the left-hand side is equivalent to the expression Jn (B - py)/o of (C3), with 
appropriate changes in notation. Thus, the authors of [9] regard samples of more 
than 10 items as adequately large in order for the use of the central limit theorem 
to give meaningful and even accurate results in most cases. The regular sample size 
in the WB method js m = 25 as stated earlier in this Appendix, and this sample size, 
therefore, appears to be adequate. 


The Chebyshev inequality rigorously assures us, as shown below, that the 50% 
confidence interval determining the probable error will never be larger than the 
Interval (p — J2 Op, pt ¥2 On) where OF (=o/ Jn) is the standard deviation of the 
mean p and o is the standard deviation of an individual item p,. This is true, however 


45 


ree ek Aare Ra etm RRS eR A erenes THREE EEE nen Rene EE Ore Bim pine Rm 
ae mete 


ho 
a 


coe ee en a 


Spies SS a ca ee es ink er es Re Ry 


small the sampie size n may be. This implies an upper bound of ¥2 0, or ¥20/yn 
for the probable error, which of course, is substantially larger than the asymptotic 
value 0.6745 o / ya, for large n, iniplied by (C5), and based on the use of the central 
limit theorem. 


The result stated is derived from the Chebyshev inequality in the following way. 
This inequality, [4, p. 135, (10)}, with the notation changed to be compatible with 
the present situation, may be stated as follows: 


: . a? 

Prob (IB ~ pol > b) < a * © (C8) 

where b is an arbitrary real number. Now, we ask what value b must have in order 

: that the stated probability shall be 1/2, in order to find 50% confidence limits. Hence, 
we put 


on = 1/2 ~ (C9) ! 
by? 
of which the solution is: . ~ 2 
: pa Feg = VEoME co) 


: so that the 50% confidence interval, as stated, is shorter than the interval (6 — Jz 2os 7" 

pt Jz 22. p Although this result is perfectly rigorous and independent of the sitions 
size n, it provides only an upper bound for the lergth of the 50% confidence inte:val 
or for the magnitude of P.E. 


For a normal population, the probable error (P.E.), is simply the half-length of . 
the $0% confidence interval. More precisely, suppote the population raean is X, 
in the notation of [10], and the standard deviation (for individual items) iso, Let | 
a random sample of size n be drawn, with sample rican X,. Then, the 50% confi- : 


( 
dence interval for that sample is ne C,), where C, =X, * 0,67450/ nh, C, =x, + F 
} 0,67450//n, or C, = X, - PE, C " x, + P.E., since the probable error “tor the 4 
4 sample mean in this Sick is 0.67450 he Thus, the length of the 50% confidence : H 
- interval is C, —C, or 2 P.E., or, the P.B. is the half-length of the confidence interval. ie 
See the derivation of (C3) above, where the sample standard deviation [Z{p, - p)*/ . P 
(n ~ 1))' (for individual items) is multiplied by 1/\/n to get the standard atte " 
Fe tion for the mean §, and then by 0.6745 to get the probable error, P.E., for the . | ‘ 
4 Bi sample mean. 
if 


Baer iF ou ela 6 Ce Ce 


et 
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We show here how this relation between the P.E. and the $0% confidence inter- 
val follows from (10, pp. 188-189], and we refer the interested reader to this stand- 
ard statistical handbook for further details. 


Let the confidence coefficient be @, in the notation of [10], and here we are 
interested in the value 8 = 1/2, By interpolating in Table IX in [10], we find that 
hg ord, a 0.6745. This well-known constant is the value Ae such that 


1 OXz x3 
over Js exp [— 1/23) dx = 1/2 


(see (C2). Then [10, p. 189], C, =X, — Ayo//n =X, — 0.67450/Vn, C, =X, + 
Ago/Vn * R, + 0.67450//n, and the 50% coufideiice interval is{C,, C,). The sign: 
fcanes of the term “confidence interval” itself is brought out in the following state- 
ment from [10]: “If it is asserted, whenever a sample is drawn, that the unknown 
population mean x lies between the limits C, and C, calculated for that particular 
sample, then the probability that the assertion is correct is 6”, with g = 1/2 here. 


In the application made in this report to the WB method, the unknown popu- 
lation mean X corresponds to p, or the true but not exactly known P,» the sample 
mean X, is p or the average of the n conditional probabilities P; for that sample and 
the population standard deviation o in each case is estimated or approximated by 
the sample standard deviation [Z(p, — B)*/(n — 1)]1/?, 


Also, the results quoted from [10] assume a normal population. The means B 
here are not exactly normally distributed, but they are asymptotically normal! for 
large sample size n as we have shown, and the normality assumption can be assumod 
to be good for n> 10. 
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In Section 3, pages 16-18, we used points (Hain? Kin)? (Mnax? Kea» Which were 
useful in simple tests, using very little computer time, for detecting cases in which the 
kill probability, P, is approximately 0 or 1, within a prescribed tolerance. For in- 
stance, if the general tolerance is 0.005, so that only two decimal digits can be assumed 
to be correct in the values of P,, then, for any case in which P, > 0.995, the true 
value is closer to 1.00 than to 0,99. 


In this appendix, we give without proofs, the equations for the coordinates of 
these points (h..., Kain)? (Ninax? Kmax 

Suppose we are given (see Figure D1) an ellipse such as those considered in this 
report, with center at (u,v) relative to the Chk coordinate system, semi-major. and 
semi-minor axes a and b, respectively and major axis making an angle ¢ with the h-axis, 
with ~90° < ¢ < 90°, 


k 


(yor a3 Kor 2) 
(h 


max? Knax ) 


Cory 2? Kort 2) 


(hain? Kinin) 


Figure D1, Values of h k 


min? “min? 


hoes’ Kmax for -90° < ¢ < 90° 


max? “max 
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There are obviously two horizontal tangents and two vertical tangents to the 
ellipse. We denote the points of tangency of the horizontal tangents, as in Figure D1, 
as (Mor 1? K yor 1) and Chior 2? Knoe 2) with Kuo 2 2 Knot 1 always, ie., (Dior 3° : 
ki or 2) 8 the point of tangency of the upper horizontal tangent. The relative magni- , 
tudes of h and h depend on whether ¢ is a first-quadrant or fourth-quadrant ; 
angle, " 


hor 1 hor 2 


The points of tangency of the vertical tangents are denoted as (h,. .. kK. 1) 
(Hy ore 2? Ket 2) with Ty ert 2 > Hyert 1? always, ie, (Hvert 2? Ket a) is the point of 
tangency of the rightmost vertical tangent (Figure D1). 

It can be shown by analytic geometry that the coordinates of these four points 
of tangency are given by the following equations. We do not give the proof here. 


(a? - b?) sin ¢ cos ¢ 


= u = 
oe Va? sin?¢ + b? cos?¢ 
= Dlb ' 
Kwon 7 V7 Va? sin? + b? cos? (D1b) : 


2 .. 2 A 
or 4 b*) sin ¢ cos ¢ (Dic) 


mee va? sin?d + b? cos?¢ 
Kora ™V + va? sin?» + b* cos*¢ > 0 (Did) 


vert 1 BU ~ Va? costo +b? sin?g (Die) 


2. ht) eid 
k ye b*) sin @ cos @ (wih 


i Va? cos*¢ + b? sin? 
hoor g “Ut Va? cos?¢@ + b? sin?g > 0 (Dig) 


2 _ w2 
‘i egal b*} sin @ cos ¢ 


vert 2 ae ay 
Va? cos?@ + b? sin*¢ 


(Dla) 


(Dth) 


It is seen from these equations that k,,, , and h,,., , are necessarily positive, since u F 
and v are assumed non-negative and the radicals + Ja? sin’¢ +b? cos’¢ and ri \ 


a? cos? + b? sin?¢ denote the positive root in each case. : ; 


We consider further the particular case shown in Figure D1, in which the el':pse i 
lies entirely in the first quadrant. The situation is not so simple if i ellipse extends 
into other quadrants. In Figure D1, we take ed hy rae Kain OK a ik The point ; 
(hain? Kmin) 89 defined has the following property. If (h, k) is any point on or inside “ ia 
the ellipse, then: a: ale | he vee 


h> bain (D2a) 


KE Kin (D2b) 


elt Reet ee ee 


Hence, because the elliptical coverage function P(h, k, o,, 9,, R) decreases with in- 
creasing |h], and also with increasing |kl, we must have 


— ae 


P(h, k, 0,5 0%) R)< Plains Kin? Gn» 22 RD (D3) 


for any point (h, k) on or inside the ellipse, with strict inequality in Sf since we 
must have strict inequality in (D2a) or (D2b) or both (the point (h ) itself not 


, min? Kinin 
; being on or inside the ellipse). 


Since (D3) holds for every point on or inside the ellipse, it must hold for the 
average value P of P(h, k, 0,» 0» R) throughout the ellipse. But the kill probability, 
P,. is equal to the average value, P, by (2-17), Hence, wo must have 


Py = P< Phin? Kin? Tr? %? RD (D4) 


But the calculation of P, by (2-12) or (2-13) in the general case, for the input values \ 
used in Table E3 in the present report, requires something of the order of 0.4 second ' 
or 400 milliseconds, whereas one call to the subroutine PKILL, for computing 
P(hain? Kmin? [n? %%» Rd» requires only about 5 or 10 milliseconds for three-digit 
or six-digit accuracy, respectively, on the CDC 6706 machine. 


Similarly, in Figure D2, we take h, = h, ko =k 


vort 2? Kmax ~~ Knop 2 and by an argument 
i similar to the one given above for (h,,,, k,,,,,)» we conclude that: 


. P, = P> P(r. Kays Sq > RD (D5) 
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: BCP! 2 UDR EE 


The calculation of P(h,..,k,.4,- 9%,» 7,» R) requires only one call to the PKILL sub- 
routine, 


The complications which are introduced when the vllipse extends into other 
quadiants than the first are sufficiently illustrated by the case shown in Figure D2. 


The relevant analogs of (D2a and b) here are 

ihi |b, 6a) : 

; od ei [a (Ds) 

| for sites points (h, k) on or inside the ellipse, and it is seen from the figure that 
_ (D6a) is not satisfied for all points (h, k) of the ellipse if we take hain mm ho , 28 be- 


fore. We must in this case take h, | = 0, siice in this case we clearly have h @ 0 for 
infinitely many points inside the ellipss. 


6) 1? Kent 1) 


k 


hor 1? “hor 1) 


I 

! | 
| 
Bs Figure D2, Ellipse in First and Second Quadrants | 
i \ 
j 
3 | 
$3 | 


The general equations for h 
which the ellipse extends, are 


min? Kmin? Umax? Kmax? regardless of the quadrants into 


hain * Max(h,,, ,» 0) (D7a) : 

1 Kays ™ max(k,, 1» 0) (D7b) 

i Hoax = Hyert 2 >0 (D7c) 

{ ” K ax = Kuor 2 >0 (D7d) 

| 

‘ with h)., ,» etc., in (L7a, b, c and d) given by (Dle, t, g and d), respectively. Equa- 1 
; Equations (D7a, b, c and d) are the general equations for all values of the angie i: 
¢. Only four values of g, namely, 0°, 45°, 90°, and --45°, ure used in the table givon q 
{ in this report, We give in tabular form below, tle values of hain’ min? h Bi es for : 


these values of 9° 


max(u--b,0) 


“max(u—a,0) max(u—p,0) 


max(u--p,0) 


kin max(y—b,0) max(v¥—p 0) max(y—a,0) maax(v—p 0) iZ 
he. uta utp utb utp 
Kiaax vtb vtp vta vto 


{ 
i 
tions (D?a and b) rofer to the algebraic maxima; thus, for example, max (—3,0) = 0, 
| 
| 
| 
| 
( 
| 
| 
i 
\ 
{ 
| 


: Here, p ™ ¥ (u? + b*)/2, It will be noted that the entries here for ¢ = 45° and 
@ = --45° are the sane, as they should be. From (Dle, b, g and d), it is seen that i 

Hert 1? Knor t* Mert 2° hor g? Which occur in (D7a, b, c and d) are the same for these 

i] 


i" two values of ¢. 

i 

i. Even though, in the program, u and v can be given arbitrary real values, it has been ; 

iN assumed in this appendix that u and v are non-negative, as is the case in Table E3 of the | 

i present report, Fora case in which u or v is negative, the best procedure is to determine 

i an equivalent case in which thuse coordinates are non-negative, as in the examples on ; | 

he pages 20--21. See also Example 2, page 58, in which u and v are interchanged. 
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APPENDIX E 


IDENTIFICATIONS FOR TABLE E3 


TABLE E3 


IDENTIFICATIONS FOR TABLE E3 


Weapon falls in uncor- 
related elliptical normal 
distribution with mean 
point of distribution at C, 
the origin of Chk rec- 
tangular coordinate sys- 
tem. Lethal radius of wea- 
pon is R, with cookie- 
cutter lethality function. 
Point target is at random 
point from uniform dis- 
tribution in ellipse. Con- 
ter of ellipse is at (u, v) in 
Chk system. Semi-major 
and semi-minor axes are 
a and b respectively. Ma- Figure El, Input Parameters o,,0,,u, Vv, a, b,¢, R 
jor axis makes angle ¢ 
with h-axis, positive coun- 
terclockwise, -90° < ¢ <¢ 90°. o, ™ 1 throughout the Table E3; 0 < o, < o,; 
a>b>O0;uz 0,v > 0;R>0. This is a direct table giving kill probability P, as 
a function of o, (* 1), o,, R, u,v, a, b, 9. 


Ranges for the variables are as given in Table El (also given on page ! 2), a, 
always has the value 1, a takes the values 0,1, 0.5, 1, 2, and ratio b/a takes the 
values 0.2, 0.5, 0.8, and 1 on every page. Angle ¢ takes the values 0°, 45°, 90°, 
—45° on every page. Theretore, the only variables whose values are given in Table El 
are R, o,, u, and v. The total number of cases in the table is 59,904, on pages 61~ 
111. Three decimal digits are given for the values of P,., but since the tolerance e, 
as discussed in Section 3, was generally 0.005, the third digit may be in error by 
several units, and the values of P, are reliable to only two decimal digits. 


a Met 
: gy ae 
- CT. agp, BAM dae 


daa 


Table El. Values of R, o,>U, V for Main Table (E3) 


Values 
of 


oy 


Values 
of 


Values of u 


Values of v 


61-64 0,2,0.5,0.8,1 0,0.5,1,1.5,2 0,0.25,0.50,1 
65-68 0,2,0.5,0.8,1 0,0.5,1,1.5,2 0,0.25,0.50,1 
69-72 0,2,0.5,0.8,1 0,0.5,1,1.5,2 0,0.25,0.50,1 
73-80 0.2,0.5,0.8,1 0,0,5,1,1.5,2,2.5,3,4,5 0,0.25,0.50,1 
81-88 0,2,0.5,0.8,1 0,0.5,1,1.5,2,2.5,3,4.5 0,0.25,0.50,! 
89-96 0.2,0.5,0.8,1 0,0.5,1,1.5,2,2.5,3,4,5 0,0.25,0.50,1 

(See Note) 

97-100 0.8,1 0,0.5,1,1.5,2,2.5,3,4,5 0,0.25,0.50,1 

101-108 0,2,0.5,0.8,1 0,0.5,1,1.5,2,2.5,3,4,5 0,0.25,0,50,1 
109-111 1 0,0.5,1,2,2.5,3,3.5,4, 0,0.25,0.50,1 


4.5,5,5.5,6 


Note: For R* 2,0, = 0.2; v= 0, 0.25, 0.50, 1.50; 
But for R ™ 2, uv, = 0.5, 0.8, 1: v = 0, 0.25, 0.50, 1 


Examples 


We give two examples of the use of the table. In Example 1 .vere is only a 
trivial scaling problem. In Example 2 we must interchange the roles of h and k and 
make corresponding changes in u and v, and interpolate with respect to v. See also 
several examples in Section 3, pages 20-21. 


Example 1. A bomber flies east along the h-axis (see Figure El) and aims a 
bomb at the origin C. The bomb falls in an uncorrelated bivariate normal distribution 
with dispersions 0, = 250 ft in range, o, ™ 125 ft in deflection. The lethal radius of 
the bomb is 500 feet. A point target is at a random point from a uniform distribution 
in the ellipse shown, with center at u = 375 ft, v = 250 ft, semi-axes, a = 125 ft, 
b = 62.5 ft, with the major axis inclined at 45° as shown. Find the single shot kill 
probability P,. 


en Se a a 


a ae re 


teens. 


fast ete 


Solution. We must normalize by dividing all distances by go, = 250. The unnor- 
malized and normalized values of the input variables are as follows: 


a 
unnormalized Faas, 
0.5 


normalized 


On page 91 of the table, entering the table with the normalized values of the variables, 
we find that P, * 0.540, Since the tolerance ¢ is 0.005, the last digit is in doubt, but 
the figures 0.54 can be regarded as reliable. 


Example 2, Suppose we are given the values o, ™160,0, 8 200, etc., as on the 
first line of Table E2. 


Table E2, Data for Example 2 


Original Data 


2. Normalized 
(divide by 200) 

3. Rearraxged with 
ot > Oy 

4. Tabulated Value 0.230 
(page 86) 


Tabulated Value 


By interpolation, P, = 0.208 for data on line 3. 


58 


We first normalize by dividing distances (all variables except ¢) by the larger 
standard deviation, 7, = 200, producing the data on line 2. Then we raust rearrange 
with o, > 0,» in order to use the table in this report, producing the values on line 3, 
The vaiues of o, and a, are interchanged, and, as a consequence, also the values of 
u and v are interchanged from the values on line 2. 


(a) (b) 


Figure E2. Normalized Data for Example 2 
(a) Line 2, (b) Line 3 of Table E2 


However, a and b are still the semi-major and semi-minor axes respectively, of the 
ellipse, so are not interchanged, but a consequence of the rearrangement is that ¢ is 
changed from 90° to 0°. A consideration of (a) and (b) of Figure E2 will convince 
the reader that these figures represent the data on lines 2 and 3 respectively, and that 
the cases are equivalent. 


However, the value v = 0.7 is not given in the table and so we must interpolate 
between the values v = 0.5 and v = 1 which are given (see lines 4 and 5 in Table E2), 
The results from the table, page 86, are P, = 0,230 for v * 0,5 (line 4) and Py, # 0.175 
for the data on line 5, with v= 1, Interpolating linearly, we find that for the data on 
line 3, P, = 0.208. Here the linear interpolation over a relatively large span, v = 0.5 to 
v= 1, introduces an additional error which is difficult to assess on theoretical grounds. 
However, this case was computed directly by the DJH program, using the original 
data (line 1, Table E2). In the computer program, o, is not required to have the 
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value | or to be not less than o,, although such restrictions were imposed in computing q 
Table E3 for this report. The direct computation, with a tolerance e = 0,005, gave a : 
; result of B= 0.210897 +, rcunding to 0.211 and thus differing by slightly less q 
1 than 3 units in the 3rd decimal digit from the result obtained by interpolating in 
| Table E3. Since the overall tolerance was ¢ = 0.005, the interpolation introduced a 
7 only a small additional error in whis case. 1 y 
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FORTRAN LISTING OF DJH PROGRAM 


EXPLANATION OF FORTRAN LISTING FOR P, 


Routine LDBSMP evaluates the integral 


oe fre Oy: 0,, R)r dr do 
with 


IhnltR ikl+/R VRT = Th 1 /x? — y? 
P(h, k,0,,0,,R) * see Y ire —~— (— +—)| dy dx 
(kl JRF — Gxt)? 2\o7 0? 


(See page 6)* 


evaluated by subroutine PKILL, where 


h*u+rcos@ cos¢ ~ crsin@ sin ¢ 
k=v+rcosé@ sing + crsin@ cosd (See (2-11)) 


Identifying the other variables 


a _ is semi-major axis of the ellipsc; a> b > 0 

b is semi-minor axis of the ellipse; a > b > 0 

ec isb/a; O<c <1 

u___is the abscissa of the center of the ellipse along h-axis 

vy __ is the ordinate of the center of the ellipse along k-axis 

is the standard deviation in the h direction; o, > 0 

o, is the standard deviation in the k direction; o, > 0 

Ris the lethal radius of the weapon; R > 0 

@ is the inclination of the major axis of the ellipse to the h-axis measured 
positive counterclockwise from h-axis, -90 < @ < 90. 


Routine LDBSMP calls subroutine ARIEL in the following manner 


CALL ARIEL (u, v, a, b, o,,°,,¢, R, Py, tolerance) 
with tolerance set to 0.005, Its basic function is to set up the input and locate the 
output, 


*See also note, page 29, regarding h and k, 
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Subroutine ARIEL calls the double integration routine DBSMP to compute Pat 
When 0, 7 F% the routine is much faster because the circular coverage subroutine 
CIRCV can be used instead of the elliptic coverage function subroutine PKILL, (see 
page 32). 


Subroutine DBSMP is called by 


CALL DBSMP (xu, xf, yu, yf, eps, N, g, PBXY), 


xu__is the upper limit of integration of x 
xf is the lower limit of integration of x 
yu isthe upper limit of integration of y 
yf is the lower limit of integration of y 
eps isa tolerance which is used to terminate the integration, when 
I 


ee I, < eps/(1 + yu — yf) (See page 23.) 


where 


I, approximates P, with partitionings of the (x, y) (or r, 6) intervals 
of integration and n, is a refined partition of n,. 


N is the maximum number of subdivisions of an interval of integration 
in either the x or y direction. 


g is the location of the value of P,. 


PBXY is a subroutine which evaluates the integrand at any point in the 
interval of integration. 


Subroutine INLINE is used by DBSMP, 


Subroutine PBXY-—3 replaces PBXY in DBSMP when o, =0,, i.6., when CIRCV 
is used instead of PKILL. 


Subroutine CIRCV is called in the following way 


CALL CIRCY (R, d, 1, P, IR), 
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where 
R 
d 


P 
IR 


is defined above 

is the distance from the origin to a point of integration (h, k), ie, 
d= Sh? +3, 

is the location where the value of P is stored 

is an error indicator which is not applicable here, 


CIRCV uses the special functions ERF and ERFD in tabulated form. They are 


defined as 


ERF(t) « —2 f ay 
vr Jo 


2 2 
ERFD = — ¢-t 
Vr 


fr 


Limitations: 


Iho, < 600 iki/o, < 600 


1/15 < 0/0, < 15 (See page 29) 
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PROGRAM LDASMP COUTPUT ) 
COMMON/ONLY1/FR(1025) ,EPSLNPTS 9N3yN5,N7,I0P5 
ToPs=14 j 
91 FORMAT ( 1HD,6E22.15 ) i 
5 97° FORMAT { 1LH1 ) 
PRINT 97 
PI=3e1415 92653 58979 ; 
UL= 140. ’ 
10 Al=2006 
AI=1006 
$xi2160. 
S¥3=200. 
P3=(90.4*PT)/180~6 
15 R3=300- 
EPS1=. 0001 | 
CALL ARIEL € Ute V19A1 983 9SX1,SY3S»P3, RIGS ,EPSL) | 
PRINT Dlg ULeV1eAL e839 SXL SVS yP3_RS96S 
CALL EXIT 
20 END 


SUBROUTINE ARTEL (UUs WV pAg 8 ySIGXySIGY) ANGLE, ARKILL)GS,EPS4 ) | 
EXTERNAL PBXY »PBXYS 
ySOHHON /TAOL/ Uy Vg PHI, CEE pSIGMAX» STGMAY »RKILL» SNPHT sCSPHI 
5 Vavy 
RKILLZARKILL 
SIGMAY2SIGY | 
SIGMAX=SIGX 4 
PI23.14159265358979 { 
10 CEESAZA 1 
TWOPT =2,*PT i 
PHI®ANGLE 
SNPHI=SIN(PH1) 
i CSPHI*C0S (PHI) 
415 XUSTWOPT 
EPREPSL#PISASA 
E4MEP { 
XL 20 
YUsA 
20 YL=0 
IF ( SIGX.EQ.SIGY) GO TO 1104 
CALL DBSMP ( XUyXL YU, YLyE 1) NMAX,GS,°BXY ) 
GS=(2e"GS)/(THOPI*A*A ) 
GO TO 1154 
a5 4404 CONTINUE 
UsU/SIGX 
VaV/SIGX ; 
| SIGMAX=4,. 
| SIGMAY21. 
| 30 RKILLERKILL/ SIGK 
Wee YUsA/STGX 
| ELSEPSL*PT*y UF YU 
| CALL DBSMP ( XUgXLyVUy YLe 1» NHAXyGS,PBXYS ) 
GSeGS/¢ PI¥C CA*AD/CSIGX*SIGX) ) ) 
bee 35 4451 CONTINUE 
| RETURN 


ENO 


a a a ee OY 
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SUBROUTINE DBSHP(XUP,X LO, YUP »Y LO,EPSR, NMAX,EVALFNy PBXY ) 
GONMON/ONLY1/FR(1025) pEPSL»NPTS 4 N3yN5yN7 yIOPS 
PIm3.1415 9265 3589 79 
NMAK=1025 
5 N304 
NBs 
N70 
EPSL®EPSR/ (1 »+(YUP=YLO)) 
: Hz (XUP=XL0) £2 


cee me oe telat nants 


| 10 Cc 
c FOR 3 LINES 
Cc Se 2 
Iy=0 7 
: Na3 a 
15 HK=a(YUP=YLO) /2. ; 
DO 100 JJet,3 : 
TYaT Ven { 


FYZYLO+FLOAT(C JJ=2) *HK 
100 CALL INLINEC FY, Judy IV) POXY, H » XLOy NMAX ) 
20 CONS 22FR (1) +FR (3) : 
OL VAL (HK (CONS24¢4 6#FR(2)))/3, 


SAVE LINES ALREADY INTEGRATED ANO MULTIPLY BY 2 


aang 


25 150 PLOwf, ; 
DO 200 Je3,N 
Lan=J+2 : 
NUMe2*L=4 ' 
FR (NUM) 2FR (LL) i 
30 200 PLOSPLOFFR (L) 
PLOs2.*PLO 


ce nmin te eee nm a a ain he ea el IP aE le De 


COMPUTE ADOITIONAL LINES ANO MULTIPLY BY & 


35 PLN#0. 
Na24Neot 
HK=HK/ 2. 
O00 300 JJm2)N,2 
Tysry+d 

40 FYaYLO+*FLOAT( UJ=1) #HK 
CALL INLINEC FY, Judy IT¥>y PBXY, H 9 XLO, NMAX ) 

300 PLN®PLNGFR (JU) ‘ 

PLN@4 o*PLN 


INTEGRATE RECTANGLE 
EVALFN2(HK* (CONS2¢PLO¢PLN)) 7/3. 


aang 


45 


TEST FOR OESIRED PRECISION i. Z 
TF CABS (EVALFN=-OLVAL)~(EPSR*ABS (EVALFN))) 500 9500,400 ‘ 
TF CABS CE VALPN=OLVAL) =¢ EPSL d 500,500,490 \ 
400 OLVALSEVALFN i 
IF ( NeLTeNMAX) GO TO 150 
500 RETURN 
55 300 CONTINUE 


TF € TOPScEQsd ) «EVALFNSEVALFNG (-25/25.) *CEVALFN=OLVAL ) 
RETURN 
ENDO 


ean aan 


50 


om 
o 


eM Rete oes 


REGISTER ALLOCATION 
& REGISTERS ASSIGNED OVER THE LOOP BEGINNING AT LINE 26 


Hoe 


et 
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SUBROUTINE INLINE( YVoyJJ»y IZ) FAXY, OHOy XLOy MMAX ) 
COMMON/ONLY1/FR(1025) ,EPSLsNPTS »9N3yN5 yN7,LOP5 


INITIALIZATION FOR FUNCTION EVALUATION AT % POINTS 


Oan 


Iv=1Z 

H=0HO 
M=s3 I 
IX=4 
40 X=XLO 
, GALL PBXY( X»y VY» I¥y IX, ANS) 

| CONS12ANS 
og TX=2 
4 X=XLO+H 

15 CALL PBXY( X»y Ys IVs IXy ANS) | 
SUM2= ANS 
1X23 | 
X=X+H 
CALL PBXY( Xy VY» IVs IX, ANS) | 
20 CONS1=CONS1+ANS | 
OFR=H*® (CONS +4 o*SUM2) /3 6 | 
| 
| 
{ 
1 
| 


SAVE PREVIOSLY COMPUTED VALUES AND MULTIPLY BY 2 
25 PFO#2.*SuM2 
COMPUTE NEW VALUES AND MULTIPLY BY & 


25 
4 SUM420. ' 1 
; 30 Mz 2"*M=4 q 
H=H/2. 1 li 


anaaqag AAO 


DO 30 J22)My2 
IXslXe1 i 
X= XLO+FLOAT(J-1)*H | 
35 GALL PBXY( XK» Yy IV» IX»y ANS) 
30 SUM&= SUM4+4ANS 
PFN 4 o*SUM4 
SUM2 = SUM2 4SUM4 


40 INTEGRATE ALONG LINE ANDO TEST FOR OESIRED PRECISION 


oaoa 


FR(JJ)=(H* (CONSL #PFNSPFO)) 73. 
i IF ( MeEQ.5 >) GO TO 50 
} Cc TF (ABS(FR( JU) “OFR) -EPSLYABS(FR(JU))) Sy Ml, 50 i 
4 45 IF CABS (FR( JJ) -OFR) -EPSL > 40,40, 50 
a 50 OFRRFR(JJ) 
IF(MeGE.MMAX) RETURN 
TF (TY eNE ot eAND«AIX cEQNPTS) L¥=4 ; 
GO TO 25 .: 
' 50 40 IFC(IY.EQ.1) NPTS#IX 4 
: N72N74¢M ; 
tl IF ( NOLT7eM ) NI2M 
; N5aN5+2 
: IF € IOPScEQeL ) FRO JI) BERE JUD 4110/1509 PCF ROIS) =OFR ) 
55 RCTURN 


| 

| { 
! END | q 

| 

! 


: aes 
rear aot ip RRO R TNE PRR A ited ie ie sh eels! 


SUBROUTINE PKILL (RRRy SAR», SYR» FHReRKR, RESULT ) 
DIMENSTON NN(6) XI (24,5 ) gHI (E496 } 


c THIS PKILL IS FOR 3 OIGIT ACCURACY. 
¢ 
5 REAL K 
INTEGER OG 


DATA WNMEG/ 29999997 3/5 TONMER/1 69999995 /9A/5 0 452 / pASTARSS 874/,AL/% EF OOL 
We 8927 A5/5 0233s pSQRT2/ 1041421 36/ »NN/ Gs 6512526) 209247, TWSQPI 41.1283 FF O02 
B792/ pSQZ02/4 0707106784, TOLS AL E077, Eros - 
10 CXT/= 6 93246951 9 066120939 9~ 62 386191992 23861919y 066120939, 93246951, KF O04 
D18*0 .0,~.9602 986» 079666648). S25552h1,—- cL BSH3HG4, 0 19363466,.5255 FFOOS 
: E3248 » 679666648 » 096028986) 1650 0 09-6 98156063) — 290421 726,~2 76990267,= EF ODS 
' Fo 58732795 _— 036 FASLS0»~ 6 125 2336 19. 12823 342 5 0 3676315050587 $1795, 2.769 EFOOT a 
GIO267 po 904247265. 98156 O63, 1270 639% 0 969400933 — 0 94H87S02,-. 86563220, FFOOR 
15 H= 07584 0b) 0 61 7E7E2Gy 0 HS O016 78, ~ 0281603555 ~ 9501 2520E= 0152980225 EFDND 
TL 0E+04 9. 281603555 6458016785 62 TE7E Zh, » 755408 U1,» BESE9120 9 94487502 EF DAO 
S96 989400 93) 840.05 ~ 5 99312 8B Dy =o 96397193 po DL 225HG Sy Mme BIGLIO9T 9 o 7HH EF OLY 
K33191 5) =. 63605368, — 231086700 p= S737 06095-2227 78585, 0 POS2ES LL E<045. EFOL2 
L 76526524 E-015 6227786565 5 6 373706099 0 521086700, « G3605368,.74639191,.83 EF OLS 
20 M9216 97 909223443, 696397193 9.99 IL 2860 40.05%. 995286722,%. 97472856,= EDLY 
Na 93827455 ~. 866 UL 6535-6 82000199) 0 7U9L 24299 — 2 64809365) 6 SKSH 2147, = EF OLS 
06433793515 ~-0315046268,— 01911188 7) e B4OS6B9STE=<01,y , 64056893 E-025.1951 EFOLG 
PLES? 9 o TL50426 8 ye GXI79S5 1 yo SUSH 2LUN7 , OH 309365 » oe 7HO01 2015970 820001995. FOL? 
QB 86425535 .93827455_ 097472856, . 99518722/ EFOLS 
25 DATA WIA 6171324495 0360 POLS 70h 6791393, 467913593, 636076157, oA 713246 KF FOLIO 
ADe 184000, LOL 228545 eo 22258103, 0 34370665, . 36268378, .SE208378, 03137 FFD 20 
80665 » 62223810 3p 010122084, 1E 700 Oy 47175 396E=015 0106939335 oA 6T07833, EFO2t 
Ce 2O0SLE7H3y 233492549 24914 705s oe ZHDLK 705 5 eo 2334925 po 20816743, 16007 EFDZ2 
D833» 610699983 , 4717SSZGE=b1, 129000» oe 2PL52459E=01 » eG 22535 PHL $029.98 EF DZS 
30 ELS 681 26=015522462697 564495 95999 1691565250 L4260342 019945064 961894 EFOZh 
FSO6Ly 616 260342 » 4269856 525 6 1495 9599 12462897 po 95 158512E=02,.622535 EFOZS 
G2GE01 5s 27452 459E 01989 O Oy oA 7ELKOD7TE ML, oh O6OLK SDE =1 9 .62672048E=1, EFOZS 
He 8327674 26915 610199012 » oe 11819453, o L316 8864s 142096115 016 917299,015 EF O27 
12753399 0615275 3395 08499729595 0142096119 L31GB8H GS yo 11819453 yo h0193012 FFOZA 
35 dp 0 832 7E7K2E~2 4 0 G267204 BE W156 40 GOLK SOE 1, AT E6LGOD7E HL, HDs Oy oA 23K12 EFOZ9 
KS 0E= 015028531 TADE “019. U427 PHSGE~ 02 y  S929SSOSEO4, 0 7334648LE=01,.86 EF OSD 
4.19016 26-014 5697621865260 15010744427 5 011 5505675 0 12167067,0425837465, EFOSL 
M42793820 5012793820) L258 3746, 0 S216 7067» oLL550867, ol OTK 27 ye D7HL 46 EFDSZ 
NB2E=015. 86190 162E-015 0733664 81 E902» 59 298585E$01, e4hO774S9E=01,020 FFOSS 
4&0 O53138GE-015.42341230E-04/ FFOSY 
RERRR 
SX2SXR 
SY*SYR 
H=#HHR 
45 K#RKR : 
H2aH*HeKeK . 
TF (SX-SY) 1,2,3 
4 SsSY i 
GO TO 4& : : 
50 2H 2 SORT(H2) : i 
KrO, ‘ 
3 SeSx ‘ 
4 TEMPSR=A5%S : 4 
IF(TEHP) 556.6 : 
55 5 TEMP=R+AL*S 


fla at 8 le Mea Se a ere 


PPS porn. 


120 


taro 


- eo 


BARR 


LF CTEMP* TENP=H*HeK4K) 858)9 
IF CTEMP* TEMP~H*¥ HeKFK) 59757 
P41. 
60 TO 75 
P=0. 
GO TO 75 
LF (R=H+AL*SX) 858510 
IF(R=-K+A4L*SY) By Spit 
SX225X*SQRT2 
65 SY2"SY*SQRT2 
- TH SY 4K*SX-A44SX 248Y 
LFCT) 13513,42 
12 IF(R*R*(SX*SX4SY#SY) T*T) 8,8,13 
13 ASX#A*SX 
70 ASY*A*SY 
TEMP £R -H “ASX 
TECTEMP) 14,154,145 
44 GOr0. 
GO TO 46 
75 15 GO=TEMP/R 
16 RARYSR<K=ASY 
IFURARY) 175917 )18 
17 GOP=0. 
60 TO 19 
80 18 GOPSRARY/R 
19 ITF CH@ASK) 20,20,24 
20 Gi=te 
GO TO 22 
24 Gix(R=HtASX) /R 
85 22 TEMP=K-ASY 
LF(TEMP) 23,523,246 
23 GiPai. 
: GO TO 25 
26 GLP=(R-TEMP)/R 
90 25 YI*#G61-G0-G1P+G0P 
IF (YI) 27,27,26 
26 RARY=SXx 
TE MP2 
T=Gi 
95 U=60 
5XeS¥ 
Hak 
Gi=G1P 
G02G60P 
| 100 SY*RARY 
Ke TEMP 
GiP=T 
| GOP2U 
| , 27 GiGOP=G1"*GOP 
105 G61 PG0*G1P*G0 
IF (G16 0P-G41PG 0) 52,32,28 
| . 28 TFC ABSCYT) = TOL) 349 34,29 
| 
| 
{ 


60 


+ Ow cd] no 


29 Ux69461460P4G1P 
T22.* (G0*G04G 146146 OP* GOP +61P* G1P) -UF0 
110 IF CT) 30,32,32 
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30 SS=64.°G1G9P*G1PG0 
IF(SS<-T*T) 31932) 32 
34 RARY=SX 
TEMP=H 
145 TaG1 
u=Go 
SXaSY 
H2K 
GisGiP 
120 6C=60P 
SY=RARY 
K= TEMP 
GiPsT 
GOPsU 
‘925 32 E41 = SQRT(G1) 
( EO = SQRT(GO) 
EN«(E4-E0) *(. 2 7*R/SX4(R/SY) 74.05% ABSCR=K) /SY4+100)) 
: SX2eSX*SQRT2 
: SY¥2eSV*SART2 
; 430 AS Ys A*SY 
IF € EN=200 ) 44940, 40 
; 40 OG24 
Go To 50 
61 IF € EN=e7 ) 43p42_42 
' 435 42 0623 
GO TO 52 
: 4&3 IF (¢ ENe.5 ) 45 44, 4u 
44 0622 
GO TO 50 
' 440 45 0624 
50 IAs. 
Test 
TF(H) 52,51)52 
54 IH=2 
£45 GO TO 53 
52 THet 
53 TF(K) 56,55954 
54 IK=4 
GO TO $6 
150 55 Ikx2 
56 Px0. 
N=NN(0G) 
00 74 ImiyN 
TIUM(EL“EO)* (XTC, OG)44.97204E 0 
} 455 VaTIJ*TIJ 
ji UsR* (1 0-V) 
i XLav* (2.<V) 
i 


Sk a ree 


ee ene 4 


2 


Ree ae eo pn 


pee ee ae 


‘ 
ert A oe te 


SS ee eee 


XQ 2 SQRT(XL) . 
XL2XQ= (ASY@K) /R - \y 

160 Q2XG- C(ASY+KI7R lj 
XI Je (HU) /SK2 { 
KeoK7 JX ae: 
Fis €&XP(X)*TWSQPI = 
GO TO (57,58) »TH 

165 57 YIde (H#U) /SK2 


XzeV¥IJ*VIJ 
F2 = EXPCX)* THSQPI 
FIJ=FitF2 
GO TO (59,72,69),T8 
170 58 FIUs2."F1 
GO TO (59y72,69),TA 
59 IFC(XL) 60,571,721 
60 Dite(KeR*XQ)/SY2 
IF (012-A*SQ202) 62,614, 64 
175 61 O112ASTAR*SQ202 
62 Pt = ERFNWL(0412) 
63 GO TO (64970) »tK 
64 DOLe(KeR*xX0)/SY2 
TF (DO1+A*SO202) 65,65.66 
189 65 0D01*-ASTAR*SQ202 
66 P2 = ERFNWL(004) 
TF(O02) 67,67,68 
67 2 = Pi + ABS (PP) 
GO TO 69 
165 68 2=Pi-p2 
69 P=Pewl (1 ,0G)*FI Ue Z* TIN 
GO TO 74 
70 P2=P1 
G0 TO 67 
198 71 PIsWNME& 
TAs2 
T1622 
72 ITF(Q) 63973,73 
73 Z=TONMES 
195 TAs3 
Tez3 
GO TO 69 
74 CONTINUE 
P#P¥R* (E1-E0) 7(SK2%4,) 
200 75 RESULT=P 
RETURN 
END 
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SUBROUTINE CIRCV (RK, OCy KR; PV, IR) 7: 
IR = 0 


IF (RK) 105 20, 30 
10 IR = 7? 


= 
i § 15 PY ® 1. : 
| RETURN q 
| 20 Pv = 0. : a 
: RETURN 

} 

| 

| 

| 

| 

! 


ho IR = 8 

1 49 Go TO 18 

100 PV = 4, , 
RETURN 


30 A ® RK 
: 15 B = ABS(DC) 
i GCs 0. : 
IF (KR) 4D, 220, 50 


WLOG = A=B 
TF (WLOC # 6010964103) 20, 20, 90 


90 IF (WLOC = 665446793) 110, 100, 100 
110 3A = A*B 


| KFLAG = 0 
| % ® SIGN(4e,WLOC) 
IF CBA = 100) 1205 £205 240 


20 


IF (A = 665446793) 222, 
222 IF (8) 290) 290, 224 
224 OFF = ASA 
i 30 SA = B*8 
: IF (COFF/SA = OFF = 40.) 230, 230, 290 
; 230 OFF = A/B*,.5 
1 A @ (1.°B) "OFF 
B = (4.+8) *OFF 
35 KFLAG x-4 
BA = A*B 
WLOG = A~B 


100, 100 


S = ERFO(.70710670*A) *ERFD( .70 71067898) 
40 DE = A*A + O*B 
DD = (2. /0E*BA) 942 
| @ = ABS(WLOC* (AeB) £08) | 
20 = ,78539816*S | 
WO = d-Q*70 
45 Z2NM2 = 20 | 


W2NMN2 = WA 
STGB = 0. 


IF CKFLAG) 130, 160, 160 


G ® Cel. 

TF (C-1000.) 140, 250, 156 

440 Z2NM2 # (.5/C*8A)4%*2 * Z2NM2 : 
OFF = DD*WENM2 


W2NM2 = (OFF = oS/C*OFF) = (He /DE*C + 16) *Z2NN2"Q 
SIGB = W2NM2+SIG8 


60 


65 


70 


7 


95 


100 


105 


410 


IF (ABS(H2NM2) = . 0000005) 1505 150, 130 
150 PV = WO + SIGB 
RETURN 


160 SIGA = 0. 
170 C = Ctie 
IF (C-1000.) 480, 200, 200 
480 Z2NM2 = (BA/C*.5)**2 * Z2NM2 
SIGA = Z2NN2+SIGA 
OFF = DD*W2NM2 
W2NH2 = (OFF - «S/C*OFF) @- (ue/DESG + 46) *Z2NM2*0 
SIG8 = W2NM2 + SIGB 
IF (AGSC(W2NM2) = 09000005) 190, 190, 170 
490 IF (ABSCZZNM2) = 20000005) 200,200,170 
200 OFF = 1.-Z0-SIGA 
PV « (CHO + SIGB)*X + OFF)*.S 
RETUGPN 


240 OFF = ABS (.70710678*HL 00) 

SA = ERFOCOFF) 

SB 2 -ERFCOFF) + Le 

S$ = & *BA 

SQ s SORT(2."BA) 

SIGA = SQ*SA 

W2NM2 = o25*SA/SQ 

Q = (A+8)/S0% S8%*. 70710 678 

sree = Q 

DE = ABS (CAB) /S* WLOC) 

DO = WLOC/S*WLOC 
250 C 2 Ctte 

IF (C-1000.) 2605 2605 2808 
260 WLOC = Q*00 

Q = DE*W2NM2 + eS/C4WLOG ~ WLOC 

SIGB 2 Q + SIGB 

Z2NM2 2 (C*2. = 1e) *W2NM2 

SIGA = Z2NM2 + SIGA 

W2NM2 = (6540) 704046) *Z2NN2/S 

IF (ABS(Q) - 20000005) 278, 270, 250 
270 IF CAGS(Z2NM2) = 20000005) 260, 260, 250 
280 PV = ((4.eSIS8)"X = SIGA/S + 1.) * 05 

RETURN 


Cc 
290 OFF = .70710678°A 


SA = ERFOCOFF) 
SB = LeERFC(OFF) 
OFF = 1.-B*8 
WLOG # SQRTCOFF) 
W2NM2 # o35SSSSS9/A*B/WLOC 
Q = SB/WLOC 
S168 # @ 
SA « SA*B/OFF 
$8 = 6*8/0FF 
OFF = SB/A*¥2 
300 C = Cri. 
IF (G-1000.) 310, 320, 320 
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| 310 WLOC = Q*S8 | s 
4 ’ Q = SA*W2NM2 + WLOC*.5/C = WLOC F 
STGB = Q+eSIGB : 
F | W2NN2 = (26*C-10) #W2NM2*OFF*(C 4.5) / (CHL) 
a4 | 115 IF (ABS(Q) = .0000005) 320, 320, 300 
320 PV = 4.-SIGB 


: RETURN 
i END 


ce km ee tate et ee 2 


2 sare Scent ee 


poet ee te 
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Elina dalek nh abi einaenh | od 1 saestmamamee * + “ON SAUE- ERR IAT Ecos aee LPs gae” ene “om 
i ire Pee bolt Soba tab EES iS DON re tre” 


jeomeomneeeenst 18s oem isthe Mat Desens ease! 90 Saree ang Se at 
anemaenees am st A aneEee eae se meee ay. —eeeTvONeY 


10 


15 


20 


35 


40 


4s 


50 


55 


FUNCTION ERF(X) 
OIMENSION TERA(415)5 TERD (451) 


KF OOANFKNFSANFP WAV FaFranr 


GaN OentOw el an m= OOInOFOUSan+ Dar OuUFranre 


DATA (TERACI), 222,125) / Ley 


13*. 99999999, 


6*. 99999998, 


5* 299999997, 9%.99999996, 2*.99999995, 2%.99999994, 


2% 699999993 5 2%.999999925 6999999919 099999990, 2%, 99999989, 


0999999865, 
» 99999962, 
299999974, 
299999962, 
699999965, 
699 99920, 
099999885, 
99999837, 
099999768, 
0939999672, 
99999540, 
099999356, 
« 99999103, 
«99998757, 
699998286, 
«9999 76467, 
«99996787, 


«99994760, 
699992931, 
959990508, 
«99987316, 
699983131, 
099977671, 
09997 0584, 
099961629, 
+ 99949665, 
099934620, 
099915478, 
0 99891242, 
0998607225 
0 99822444, 
099774745, 
«99715485, 
099642346, 
© 99552488, 
© 99442628, 


DATA (TERA (T) 


©9921 5622, 
699034681, 
«986647420, 
036537850, 
098249279, 
097 864284, 
097454704, 
¢ 96951621, 
096368407, 
e 93668725, 
294901604, 
» 94001503, 


© 99999987, 
© 99999961, 
699999972, 
099999959, 
099999940, 


099999924, 


099999877, 


699999825, 


DATA (TERACT), 


«99999751, 
« 99999649, 
«99999507, 
¢ 99999311, 
© 99999042, 
99998674, 
09999817 3, 
© 99997495, 
099996582, 


1"126,223) / 99995632, .99995358, 99995067, 


099994434, 
699992699, 
099989938, 
© 99986567, 
099982454, 
© 99976397, 
999966934, 


© 99959305, 


099946942, 
699931149, 
099911073, 
¥ 99868682, 
099853728, 
998613715, 
099763861, 
299702053, 
-99625410, 
¢99532227, 
099417953, 


TH 220,321) £.699308990, 


099162075, 
298994516, 
298769094, 
098500283, 
e9B18LUGb, 
© 97603609, 
097360263, 
096844350, 
96237290, 
695537616, 
094731240, 
93806516, 


2 99999966, 
« 99999979, 
99999969, 
« 993999956, 
¢ 99999936, 
« 99999908, 
0 999996868, 
099999812, 
099999734, 
«99999624, 
299999673, 
© 99999264, 
«99998977, 
«99998586, 
099998053, 
0999973335 
«99996365, 


299994090, 
099992042, 
© 99989335, 
99985776, 
099981118, 
099975054, 
099967198, 
099957071, 
099944082, 
«99927506, 


©9990 6437,. 


09987 9661, 
699846423, 
699804594, 
699752529, 
09968 8046, 
6 9960 65615 
099521241, 
099392257, 


0 99147249, 
« 9895 2454, 
098719028, 
© 964460704, 
© 96110492, 
097720664, 
097262812, 
09672 76758, 
696409351, 
0 BHIGS244, 
094506144, 
093606312, 


099999985, 
99999977, 
99999967, 
99999952, 
699999931, 
©99999903, 
99999858, 
99999796, 
99999715, 
99999598, 
0999994 36, 
099999214, 
09999890 a, 
099998492, 
«99997926, 
99997162, 
4999961 34, 


099993725, 
© 99991559, 
299968695, 
°99984941, 
699980029, 
099973640, 
699968371, 
299954724, 
»99941060, 
699923686, 
299901620, 
099873766, 
099638783, 
© 99795063, 
099740702, 
0 99673KH1, 
099590638, 
2 29489200, 
099365565, 


299211103, 
«98909050, 
«98GETIB7, 
«98379046, 
+980 37559, 
097634838, 
097162273, 
96610515, 
«95969503, 
95228512, 
094376220, 
293400794, 


699279043 


099999983, 


© 99999976, 


99999964, 


099999948, 


© 99999926, 
© 999998935 
6 99999808, 
099999784, 
2 99999694, 
«99999570, 
«99999397, 
«39999160, 
099998835, 
«99998392, 
099997791, 
0999969480, 
699995890 / 


0999933839, 
699991048, 
299988026, 
099984060, 
099978681, 
699972151, 
¢ 99963450 » 
099952257, 
099937928, 
099919679, 
« 99896551 5 
© 998673 B75 
099830795, 
0997851415 
099728361» 
099658215, 
0 99571905, 
099466372, 
099337823 / 


9 099247932, 


699073598, 
098864058, 
© 96613459, 
198315259, 
0 97 9621785 
097546202, 
0 97058869, 
«96489786, 
© 95429657, 
«95067330, 
094191372, 
0 ISLO9AGS » 


60 


65 


70 


78 


60 


«92973419, 
© 91805040, 
290483743, 
«68997067, 
67332616, 
2 8547 8621, 
083423150, 


OATA (TERACT), 


379690821, 
077066806, 
074210096, 
07411415563, 
06778601195 
064202933, 
«60385609, 
256332337, 
052049968, 
047548172, 
042839236, 
037938205, 
«32862676, 
227632639, 
022270259, 
016799597, 
0441246292, 


92751363, 
91553388, 
090200040, 
«88678789, 
086977330, 
«68083402, 
82987029, 


Te322,615) / 181156356, 


092523554, 
©91295551, 
«89909620, 
¢ 88353301, 
«86614359, 
«84681050, 
0 82542365, 


«92290013, 
091031398, 
089612364, 
086020507, 
86243611, 
«8427007 9, 
0820890812, 


092050518, 
«90760829, 
289304233, 
067680310, 
2 85 864995 » 
«63650807, 
261627102 7 


080676772, 80188283, 


«79484325, «76668732, «78463985, «77610027, 


676514271, 
073.6103465, 
070467806, 
0670864006; 
063458583, 
«59593650, 
055493925, 
051166826, 
+46622512, 
241873870, 
636936453, 
31826350, 
26570006, 
024183989, 
015698703, 
210126059, 


0759523765 
o73001043, 
0698610394, 
¢ 66378220, 
eb2704644, 
238792290, 
054646410, 
«80274967, 
«45688669, 
04N90 0948, 
035927665, 
0307680075 
e 25502240, 
o2n093S as, 
eo 14586711, 


«75381075, 
e72362161,5 
069143312, 
«65662770, 
061941146, 
057984581, 
«53789863, 
«4 9374505, 
«44746762, 
039920598, 
034912599, 
029741822, 
024429691, 
018999246, 
el Sh758355 


074800328, 
071733678, 
068466555, 
064937669, 
0611681275 
o5 7161576, 
0529263624 
048465539, 
243796909, 
«38932970» 
0338908155 
028689972 4 
« 23352192, 
017902181, 
012362290, 


0090078126, .078857720, «067621594, 


0056371978, 0045111106, 4033841222, 022564875, 1011283416, 0./ 
DATA (TERD(T)» Imd,y 78) / .1821S3US9E~08, 1981 6RG2E-08, 


Rn Om ee ek 


90 


100 


108 


410 


Peer OR aero. Oe Oe ea 


SONOS NW KOMBAT AVE CSSEnrKVeunre BSIOnNPAREUNe Van rVvseanr mOeaverus 


0216765 96E- 08, 
03097 0439E=-08, 
o441076K7E-08, 
06261677 2E= 08, 
2 6460 8977E=08, 
oi 2498993E-07, 
e757 RGBKE=07, 
o246320K1E 07, 
0 S441 S471E-07, 
047925547E= 07, 
0 66529674E~07, 
09206 06948-07, 
0 12696235E-06, 
oA 74591 35E~ 06, 
o23928327E~ 06, 
0 32669796E~06, 
04451 6637E-06, 
© 60428629E~ 06, 
© $81766120E= 06, 


oL2792741E= 05, 
o17171961E~ 05, 
0 2297 6636E - 05, 
 30645257E~ 05, 
e4OPK27H9E “05, 
© 33994268E= 05, 
743272116 =05, 
093923243E= 08, 


e237 06410690 by 
0338430 33E-08, 
4816021 0E-0 8, 
e68315260E=0 8) 
+965 9559BE=0 85 
13614673E-07, 
19427901E=07, 
26767 BK1E <0 75 
e373 9SUAGEM0 7, 
+52036639E=0 7, 
27217 B9Z3E=075 
©99797993E=075 
11375445 8E=06, 
1889626 0E=06, 
02587 7218E-06, 
e35324013E=06, 
©48065419E=06, 
65193709E-06, 
06814321 E~06, 


0137 73933E-05, 
1 8G7K250E~ 05, 
e24699S74KE-08, 
¢32916626E-05, 
oh3727530E "05, 
o37903503E-05, 
0764302995 -055 
L0056235E-04, 
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e25920474E-08, 
036974673E-08, 
«52874603608, 
«745% 7H38E°08, 
«105281 02E=07, 
01482 B97KE-07, 
e 2081 4HG3E=07, 
02912 G498E=0%'5 
04062 7618E-07¢ 
©56489121E+97, 
e78292207E-07, 
eLOSLG6394E-06, 
oL4S9S557E-06, 
020 4m 7846F+06, 
©2797 9245E-06, 
» 381628667E-06, 
o518686725E-06, 
» 703204738 =06, 
09499867 9E-06, 


01462 7416E-05, 
019872 328E-05, 
» 2654 5968E~-05, 
0353492 7SE°NS, 
©4692 15869E-05, 
06206 3353F-05, 
e BL BBL B9RE~05, 
010 764921E-04, 


o 2833600 2E-058, 
04038601 6E-08, 
aB73G2INU E08, 
081266 442E=0 6, 
eLL47244 SFHN7, 
ef 61 43994E-07, 
© 2264520 4E“07, 
o31662977E=-07, 
e44L3036KE-07, 
eOL3L031 SE-07, 
0869062816 -07, 
4472077 66-06, 
©16128098E-06, 
22421787806, 
oS024S971E -06, 
04122162 bE -06, 
© 56000632E~-06, 
o TEB35232ZE=06, 
010236666E-05 / 


DATA (TERD(T)» 1279,256) / cAL0Z84G5E-05, .11879066F-05, 


01595 8261E-05, 
02136976 2E-05, 
2652491 4E=-95, 
03795412 3E<05, 
oSN336887E-05, 
266351620605, 
¢67704910F -05, 
o11S821266F-04, 


ee ce ee ee 


11S 


120 


425 


135 


440 


‘4 
4 
; 


| 145 


150 


155 


160 


Ln a ee oer NO 


£65 


Paeanvneawe ah wo 


ANH HMOOnNFVUFSANr BAN Furlan r PB OMNEVSAUN PrP OAAaNOVS UN 


0 12328243E- 046, 
2 16130192E-Oh4, 
e2L037210E= D4, 
027 349351£= 04, 
03544 1831E= 04, 
0 45782082E - Oh, 
e5O9S01S7E= 04, 
» 7966.3267E=04, 
© 9680 4434E= 04, 
e12345A98E- 03, 
0 19694459E-03, 


©2236 0761E= 03, 
© 20199597E=03, 
eo SSHUQHKGE=03, 
0 44420794E- 03, 
o 5848 U7 22E= 03, 
©69082932E-03, 
v85736992E~-03, 
eL0607090E-02, 
©1308 0500E~02, 
160791 376-02, 
1497020 48E-02, 
0 24064136E-02, 
029298098E=02, 
oe S5S564B7E M02, 
e430L 3869E~02, 
o31669067E- 02, 
1 G2347K40E-02, 
0747031 76E=02, 
s892Z215S51LE-02, 


011576041E~02, 
eL 372 5650E~ 01, 
©26198805E~ 01, 
o19070402E+01, 
e2237932Z4E-01, 
0 26178475E-08, 
305247 40E=niy 
e 384780 77E=01, 
© 44105318E-04;, 
o47471879E= 01, 
054649361E-01, 
0S2711060E~014, 
o 7A732ZOHDE-01, 
061786871E-01, 
09298 7046E-01, 
oL0534307E+00, 
011893029E+00, 
213387948E+ 00, 
01850228S0E+ 00, 


01774941 3E + 00, 
619757660E +00, 
02192 3353E +00, 


01348912 7E=-0 4, 
01724276 BE~04, 
022470263" -04, 
©2918902S5F-04, 
03779560 4E-04, 


ob 870353 2E=0 by 


©62764699E-04, 
eBO494817E“04, 
©10290362E-03, 
eA TLL IOK7E-03, 
62665661 9F-03) 


0237031446 -03, 
02986859860 Sy 
«375175086E=0 3, 
¢46974632E <0 3, 
058627725E~-03, 
©72937850E=-03, 


0904 50949E-0 3, 


41184075E-0 2) 
08377730 48-02) 
«169221 366-02, 
e207 16409E~0 2, 
¢25285285- -025 
e30760230E-02, 
¢3730109ZE=0 2, 
0450 88292E-0 2, 
954327063 -02, 
65249775 E=0 25 
o TOL1B816KE=“02, 
093225623E-02, 


el 2Z081043E"04, 
eAGSO25H5F—0 1, 
11687844 6F “Ot, 
/A9VBSUGIGE=DL,y 
023280999E~0 4, 
e272 L442 E-0 4, 
0317 038S31E-0 4, 
036819865E-01, 
o4h262485 4UE~O 1, 
ehDLB740LE“0145 
0565 76979E-04,5 
64873390E<01, 
07416610 3E-04, 
086473470E-04, 
299931 799E=01, 
010859632: 40 0» 
0122540006400, 
01378322 7E+00, 
615483 7616400, 


« 14107306E-0%, 
2 18426397E= 04, 
0233961 35E~04, 
o S1146217E- 04, 
o40297636E “04, 
e51971360E= 04, 
«66 6126 74F 0b, 
«55617 765E— 04, 
010.93 6521E= 03, 
0139253 05E-03, 
el TE74231E-03, 


e 25121089E-03, 
o SIG30GSSE=03, 
o S9BIE2Z7HE= O35 
e49665360E-03, 
eGiD3BS7THE03, 
» 76992476E=03, 
0 9540277 RE=03>5 
oe L17B37T6GE=02, 
eA KSORS25E - 02, 
eh 78057 71E “025 
0217828K2E=02, 
2 265630895 ~ 02, 
o S228 O671LE 025 
¢39123473E=02, 
eo S725330GE~ 025 
« 568901472E-02, 
06827 3562E-02, 
© 81672930E=-02, 
¢ 9738 9920€-02, 


01260 8554F-01, 
01491 1571E<01, 
oi17S583088E-01, 
« 20666965E~01; 
0 24244150E-01, 
o28279S10E=-01, 
032921 681E~01, 
o 38203897E-01, 
oAhb191723E 901, 
eo 509S47T26E-01, 
eS8S6S016E-01, 
©67096878E~ 01, 
e 766260 82E“01, 
o 87229089E-08, 
09698 1951E-01, 
2 111959546E +00, 
0126234 02E +00, 
044618673 74E+00, 
¢138946171E+00, 


»1S066387E-04, 
21969161 3E-94, 
025620 500F =O 4, 
03 3227997E~04, 
4298670 7E 06, 


1S5356429E-04, — 


e7TL 107 499E 90 4, 
09106854 2E-04, 
042620929E-03, 
01478491 9E-03, 
018750262E-03 / 


DATA (TERO(L)» IsiS%o234) / .19887624E=03, «210901 B4E-03,y 


22661853 SF-03, 
« 33488688F-03, 
041997 400E-03, 
o5249971 3E-03, 


eB S41 BO7LE“03, 
981 256247E-03, 


010060558E-02, 
e1 2416485E-02, 
olLS27S07BEHN2, 
eL B73L600E~02, 
e22897S83E 0 2, 
027899886E=02, 
oS3IBRG7TOOE-02, 
o4L026681E-02, 
eh D95L2374E-02, 
oS9562287E “02, 
e7L423190E-02, 
085372378E-02, 
010271 986E-04 / 


DATA (TERD(I)» [e255,312) / «L0622108E-04, 11 0899S 0E-01, 


ef 3450 207E “01, 
eLS543K22E-01, 
ef6343482E 903, 
e21508270F 01, 
e251 79685E~04, 
e29383624E-01, 
eS4179892Z2E-01, 
039632 026E-01, 
e4SON7 O27E-01, 
oS2774996E-045 
060608644E-01, 
e69382697E Ci, 
eo 79LTZATIEW“O2, 
0 90056524E-01, 
0102118666€+00, 
011540 363E+00, 
eL300LS40E +00, 
0146006516400, 
016 343862E¢00 / 


QATA (TEROCI)» 126323,390) / .16802916E+00, .17271408E+ 00, 


eAG236999EF00, eLB73423ZE+00, 019241273400, 


© 20284406E +00, 
0224 89587E+0 0, 


129 


¢ 2062 0799E +080, 
¢ 230658336400, 


Tae seembah ia Unbenhensinashd nateaedieeiminedaiaaia tenant 


021367101400, 
e236S2012E 400, 


ath et 


Pps ee NWT On eee ek 


170 


175 


185 


190 


195 


205 


210 


215 


tate, A airshow 


OMSAN EP OBVEMNES UNS HP OSVFRCSGUN* OAONTUNE 


¢24248443E+ 00, 
026734435E400, 
2293812246+00, 
o32186671E +08, 
©3514 7462E +00, 
2382579906 +00, 
© 41510750E+00, 
o44896267E 400, 
64 8602546E400, 
o5'2016058E+00, 
0597207 40E +00, 
©59498579E+00, 
063329574E+ 00, 
067191881E+ 00, 
071061978E4+00, 
074924872E +00, 


©80604334E +005 
+84297518E+00, 
08787 A258E +004 
oQLSIBULL E400, 
© 94590063E +00, 
9765698400, 
oL0052957E 401, 
o1032Z609E4 OL, 
6 10546218E4 01, 
0 107S0662E+01, 
0109240 56E+ 01, 
0110647036402, 
01247152664 01, 
e1L2U3ZUTE +015 
eA L279279E + Ody 


ERF @ 1, 
ARGNX © be 16 
Go TO 40 
ENTRY ERFOD 
ERF = 0. 
ARGMX @ 4<8 
10 V 2 ABS(X) 
SIGN = 1. 
TFCReLT.00) SIGN = “4, 
TF (CARGMX oLEe V) 
L * 200.49V + 5 
Xi #2 ,O4*FLOATIL) 
OX = V=xXt 
dos &84-L 
IF (ERF) 20, 30, 
20 ERF =(TERA( S236) + TERDO(J) 4OX* (10-0X4X1)) SSIGN 


RE TURN 
30 ERF sCTERD(I) * Che = OX*( DK. + 2o*X19( 2 e-DK9KL) ) 


RETURN 
END 


22485483 6E+0 0, 
e273 B108K4E+00, 
© 30067728E+00, 
032912497E+0 0, 
0 35911315E+00, 
03905 86164E+00, 
o42345088E+0 0, 
©48761S525E+00, 
ob9296567E +00, 
52934346 8E+0 0, 
© 536659089:+00, 
26045216 7E+0 0, 
06429310 7E+00, 
068159626" +0 0, 
o 7202781.9E+00, 
«7587235 8E+0 0, 


0615 36635E+0 0, 
088204344 ' 40 C, 
o 88752573 E+0 by 
092453204 +00, 
09537842 7E+00, 
0984011438400, 
240219832E+01, 
01037 3633F401, 
eL06O0L41E+O1, 
« 1079696926401, 
21096235 7E+01, 
0 11094698E+0 2, 
01119276 26+ 0 1, 
ofl 2556275401, 
0112 62663E+01, 


RETURN 
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»25471302E400, 
2280377676400, 
e30764223E +00, 
033647960E +00, 
0 S6684430E +00, 
0 39667140E¢ 00, 
0431867557E4+00, 
°46635081E+00, 
o501964657E+08, 
eSSOSAR079E¢ 00, 
©37602720E +00, 
062408756E+00, 
06825 6247E+00, 
069127486E400, 
©72992188E+00, 
0766267146400, 


8246 3224E¢ 0D, 
oA6L037 0364005 
«69617662E +00, 
°92977025E+00, 
FG 1SUL SUE +00, 
© 99122100E400, 
010165522€401, 
o210432932E +015 
010682209E +01, 
0108461 308E401, 
010998593E4¢01, 
0111224696404, 
011211 606E 401, 
0112657526404, 


e26097837E400, 
0 28704457E+00, 
¢ 3147055 7E+00, 
0 34392976E+00, 
© 37466696E 400, 
4066471 3E+00, 
«4403797 9E +00, 
04 7515331E +00, 
0511034716400, 
25478697 2E 400, 
« 58548322E+00, 
262368006E #00, 
e66226628E 400, 
0700950676400, 
¢ 73954676E +00, 


oT777751L 86400 7 
DATA (TERO(T) » 18391,451) 4 a7 B7243K3E+00, .79666759E 400» 


oe 83383665E+00, 
« 869951556 +00, 
©90473087E 400, 
o9376945KE 400, 
¢ 96916756E+00, 
oe 9982A371F 400, 
010249694E+01, 
01049047 1E+01, 
eL 0702393E401, 
108637126401, 
o110327H1E 404, 
01114508 1E401, 
o11228636E 401, 
of12736616 401, 


eli 28S792E+01 / 


))*SIGN 


SUBROUTINE PBXY (XyVpIXeIY,ANS) 
COMMON /T104/U5V» PHI» CEE) SIGM AX,SIGMAY,RKILL»SNPH I,GSPHT 
Rey 
Ti=COS¢X) 
5 T3aSIN(X) 
A=U+R*T1*CSPHI-CEE*R*T34SNPHI 
BaV+R TA*SNPHI+CEE*R* T3*CSPHI 
A= ABS (A) 
Bw Aas (A) 
40 GALL FKILL € RKILL»SIGMAXsSIGMAY) ApSy ANS ) 
ANSaR* ANS 
’ RETURN 
END 


SUBROUTINE PBXYS(X pV IXpIY¥y ANS) | 
COMMON /T101/U) Vy PHT, CEE, SIGMAXy SIGMAY yRKILLy SNPH I yCSPHT 
91 FORMAT ( 4H0,SH IR=,I7y 29H CIRCV HAS BAD INPUT VALUES ) ) 
RsY | 
5 Tiscos( x) 
T3¥SINCX) 
AuU+R*T1*CSPH I-GEE*R*T3*SNPHT 
Gav +R*T4* SNPHI*CEE*R* T3*C SPHE 


| 
] 
O1=SQRT ( A*ACB*B ) 

10 CALL CIAGY ( RKILLsOLy4,ANSeIR ) ! 
| 


IF ( IReNE}O ) PRINT 912,IR 
ANS #ANS4R 

RETURN 

END 


ES ES, 


| 
| : : 
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